%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                 %
%                            ROOT FILE                            %
%                                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Run LaTeX or pdfLaTeX on this file to produce your thesis.
%  To produce the abstract title page followed by the abstract,
%  see the file abstitle-phd.tex or abstitle-mas.tex.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass[10pt]{article}
 \topmargin -1.5cm        % read Lamport p.163
 \oddsidemargin -0.04cm   % read Lamport p.163
 \evensidemargin -0.04cm  % same as oddsidemargin but for left-hand pages
 \textwidth 16.59cm
 \textheight 21.94cm
 %\pagestyle{empty}       % Uncomment if don't want page numbers
 \parskip 7.2pt           % sets spacing between paragraphs
 %\renewcommand{\baselinestretch}{1.5} 	% Uncomment for 1.5 spacing between lines
 \parindent 0pt		  % sets leading space for paragraphs

% Use the first command below if you want captions over 1 line indented. A side
% effect of this is to remove the use of bold for captions (thesis default).
% To restore bold, also include the second line below.
\usepackage[hang]{caption}      % to indent subsequent lines of captions
\usepackage{cite} % Make references as [1-4], not [1,2,3,4]
\usepackage{url}  % Formatting web addresses
\usepackage{ifthen}  % Conditional
\usepackage{multicol}   %Columns
\usepackage[utf8]{inputenc} %unicode support
\usepackage{array}
\usepackage{float}
\usepackage{graphicx}
\usepackage{algorithmic}
\usepackage{algorithm}
\usepackage{amsthm}
\usepackage[cmex10]{amsmath}
\usepackage{float}
\usepackage{amsfonts}
\usepackage{multirow}
\usepackage{rotating}
\usepackage{longtable}
\usepackage{afterpage}
\usepackage{stmaryrd}
\usepackage{eqparbox}
\usepackage{color}
\usepackage[mathscr]{eucal} % provides \mathscr (Euler script)
\usepackage{subfig}
\usepackage{calc}
%\usepackage{ragged2e}
\urlstyle{rm}
\usepackage{authblk}

%\include{rpichap2}

\floatstyle{ruled}
\newfloat{program}{thp}{lop}
\floatname{program}{Algorithm}

%\renewcommand\Authands{ and }

%\newtheorem{algorithm}{Algorithm}[subsection]

\newtheorem{mydef}{Definition}

\newcommand{\tablefont}[1]{\fontsize{#1}{#1}\selectfont}
\newcommand{\INDSTATE}[1][1]{\STATE\hspace{#1\algorithmicindent}}

\renewcommand{\captionfont}{\bfseries} % bold caption (needed with caption
                                       % package to restore boldface.)
\renewcommand{\algorithmiccomment}[1]{\hfill\eqparbox{COMMENT}{// #1}}
%\includeonly{rpichap1}  % use \includeonly to process only
                         % the file(s) listed inside the braces
\def\NoNumber#1{{\def\alglinenumber##1{}\State #1}\addtocounter{ALG@line}{-1}}

\setcounter{secnumdepth}{5}

\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}

%    \makeatletter
%    \renewcommand\paragraph{\@startsection{paragraph}{4}{\z@}%
%                                         {-3.25ex\@plus -1ex \@minus -.2ex}%
%                                         {0.0001pt \@plus .2ex}%
%                                         {\normalfont\normalsize\bfseries}}
%    \renewcommand\subparagraph{\@startsection{subparagraph}{5}{\z@}%
%                                         {-3.25ex\@plus -1ex \@minus -.2ex}%
%                                         {0.0001pt \@plus .2ex}%
%                                         {\normalfont\normalsize\bfseries}}
%    \makeatother

\date{}

\begin{document}

%\title{Distance and time matter: Host-pathogen association analysis of tuberculosis patients via Unified Biclustering Framework (UBF)}
\title{Host-pathogen association analysis of tuberculosis patients via Unified Biclustering Framework}

%\author{Cagri~Ozcaglar \and B\"{u}lent~Yener \and Kristin~P.~Bennett
%\thanks{1 Computer Science Department, Rensselaer Polytechnic Institute, Troy, NY USA}
%\thanks{2 Mathematical Sciences Department, Rensselaer Polytechnic Institute, Troy, NY USA}
%}

\author[1]{Cagri Ozcaglar\thanks{ozcagc2@cs.rpi.edu}}
\author[1]{B\"{u}lent~Yener\thanks{yener@cs.rpi.edu}}
\author[1,2]{Kristin~P.~Bennett\thanks{bennek@rpi.edu}}
\affil[1]{Computer Science Department, Rensselaer Polytechnic Institute}
\affil[2]{Mathematical Sciences Department, Rensselaer Polytechnic Institute}


%\author{
%    Cagri Ozcaglar\\
%    Computer Science\\
%    Rensselaer Polytechnic Institute\\
%    ozcagc2@cs.rpi.edu
%  \and
%    B\"{u}lent~Yener\\
%    Computer Science\\
%    Rensselaer Polytechnic Institute\\
%    yener@cs.rpi.edu
%  \and
%    Kristin~P.~Bennett\\
%    Computer Science\\
%    Rensselaer Polytechnic Institute\\
%    bennek@rpi.edu
%}

%\thanks{C. Ozcaglar and B. Yener are with the Computer Science
%Department, Rensselaer Polytechnic Institute, Troy, NY USA (e-mail: ozcagc2@cs.rpi.edu; yener@cs.rpi.edu).}
%\thanks{K. P. Bennett is with the Computer Science Department and Mathematical
%Sciences Department, Rensselaer Polytechnic Institute, Troy, NY USA (e-mail: bennek@rpi.edu).}
%}

\maketitle

\section{Introduction}
Tuberculosis (TB) is an airborne disease which is a leading cause of death worldwide. According to World Health Organization, one third of the human population is infected either latently or actively with TB \cite{WHO}. \emph{Mycobacterium tuberculosis} complex (MTBC) is the set of species which causes TB. MTBC isolates from TB patients are genotyped using multiple biomarkers for tracking TB transmission, TB control, and examining host-pathogen relationships.

Earlier studies have found associations between TB patients and the MTBC strains which infected them. Hirsh et al. showed that a TB patient's place of birth can be used to predict the geographic origin of the MTBC isolate \cite{HirshStableAssociationPNAS}. Gagneux et al. defined the population structure of MTBC strains using six phylogeographic lineages and showed that these lineages are adapted to particular human populations defined by place of birth or risk factor \cite{VariableHostPathogen}. Visual inspection via host-pathogen maps enable making inferences from patient data and strain lineages \cite{HostPathogenMaps_BIBM_CAME}. Although names of phylogeographic lineages imply an association between MTBC isolates and patients' place of birth, none of these studies combine genetic proximity of MTBC strains and spatial proximity of TB patients together. In this study, in addition to the distribution of MTBC isolates to their host's country of birth, we add genetic proximity, spatial proximity and time into domain knowledge of host-pathogen association analysis.

Multiple sources of information can be incorporated into data analysis via data fusion \cite{L2normMoreau}. Recently, there has been considerable work on genomic data fusion \cite{LanckrietStatFrameworkGDF, GenePriorGDFmoreau}. In the TB context, Ozcaglar et al. built the tensor clustering framework (TCF) to cluster MTBC strains using multiple biomarkers simultaneously through genomic data fusion \cite{TCF_MBT_BMC}. Genomic and phenomic data sources are also combined in earlier studies \cite{GenomeInteractomeLageMoreau} via genome-phenome data fusion. However, there is no significant work on genome-phenome interactions of MTBC isolates and TB patients.

In this study, we present host-pathogen associations of tuberculosis by incorporating genetic proximity of MTBC strains, spatial proximity of TB patients, and time into domain knowledge via Unified Biclustering Framework (UBF). We simultaneously factorize multiple sources of information in various forms and obtain biclusters which represent host-pathogen pairs, while keeping pathogens genetically close in order to estimate most likely mutation events, and keeping hosts spatially close in order to estimate most likely transmission events. Based on factor matrices of hosts and pathogens, we generate the feature pattern similarity matrix of host-pathogen pairs, and find density-invariant biclusters. Finally, we select statistically significant biclusters among them and find most stable host-pathogen associations. We also find host-pathogen associations within each major lineage. We evaluate biological relevance of statistically significant biclusters, confirm known host-pathogen associations, and propose new ones.

\section{Background}
In order to find relationships between MTBC isolates and TB patients, we uniquely identified them by their characteristics. We represented MTBC strains with a commonly used biomarker, spoligoytpe, and represented each patient with their country of birth. Finally, we stated the host-pathogen association analysis as a biclustering problem. Next, we give a brief background on spoligotyping, biclustering, and explain host-pathogen association analysis as a biclustering problem.

\subsection{Spoligotyping}
Spoligotyping is a DNA fingerprinting method of MTBC which exploits the polymorphism in the DR region consisting of 36 bp of direct repeats separated by 36 to 41 bp of spacers \cite{KamerbeekSpoligotyping}. A spoligotype consists of 43 spacers, and it is represented as a 43-bit binary vector, where zeros represent absence of spacers and ones represent presence of spacers. Mutations in the DR region can result in loss of spacers, but not gain. This rule of irreversible mutation of spoligotypes is also known as contiguous deletion assumption \cite{InsightsIntoBIBM2011, InferredSpolIEEETNB2012}.

\subsection{Biclustering}
Biclustering is a class of clustering algorithms which perform simultaneous clustering rows and columns of a matrix. The term was first coined by Cheng and Church for gene expression data analysis \cite{ChengChurch}. Following them, many biclustering algorithms motivated mostly by bioinformatics applications are developed. These biclustering algorithms include spectral biclustering algorithm by Dhillon et al. \cite{DhillonSpectralBiclustering} and Kluger et al. \cite{KlugerSpectral}, Statistical-Algorithmic Method for Bicluster Analysis (SAMBA) by Tanay et al. \cite{SAMBAtanay}, Coupled Two-Way Clustering (CTWC) by Getz et al. \cite{CTWCgetz}, Binary Inclusion-Maximal biclustering algorithm (BiMax) by Prelic et al. \cite{BiMax}, and Densely-connected biclustering (DECOB) by Colak et al. \cite{DECOBrecep}. A great survey by Madeira et al. details biclustering and existing biclustering algorithms for biological data analysis \cite{MadeiraSurvey}.

\subsection{Host-pathogen association analysis: a biclustering problem}
Biclustering was initially motivated by gene expression data analysis in order to group genes into subsets of genes which are coexpressed under certain subsets of conditions. This is equivalent to finding submatrices in a gene expression matrix such that the submatrix entries follows a cohesive pattern under investigation. In the TB context, the genes of microarray data maps to spoligotypes of MTBC strains, and the conditions of microarray data maps to country of birth of TB patients. The resulting host-pathogen matrix of tuberculosis expresses the association level of a spoligotype to a country.

In the case where the original host-pathogen matrix is extended or concatenated with other matrices via data fusion, we use feature patterns for spoligotypes and countries. We first extract feature patterns for each spoligotype of MTBC strains and for each country of birth for TB patients. The association level of a spoligotype and a country is calculated as the similarity of their feature pattern vectors. This final form of host-pathogen matrix of tuberculosis expresses association level of host-pathogen pairs, and is in the correct form to be analyzed via biclustering.

In the next section, we present the methods used for host-pathogen association analysis. We first give details about the patient dataset. Then, we present calculation of genetic proximity matrix and spatial proximity matrix used in data fusion. Finally, we present the steps of Unified Biclustering Framework (UBF).

\section{Methods}

\subsection{The dataset}
The NYC dataset consists of 4876 patients in the United States detected from 2001 to 2007. The spoligotype of MTBC strains and their host's country of birth are available in the dataset, along with the date of diagnosis. There are 858 unique spoligotypes in the original dataset. MTBC strains are labeled by major lineages based on their spoligotypes using Conformal Bayesian Network (CBN) model \cite{CBN}, and by KBBN sublineages using the Knowledge-based Bayesian Network (KBBN) model \cite{KBBN}. We refer to spoligotypes using shared type numbers, or SIT numbers using SITVITWEB database \cite{SITVITWEB}. If spoligotype is not assigned an ST number by SITVITWEB, then we assign a unique UST number, where U denotes unknown ST. We first filter this data such that there are at least 2 patients from a country, and at least two patients infected with a strain. After filtering the dataset, there remains 4301 patients, 311 spoligotypes, and 104 countries. Using this filtered dataset, we construct the host-pathogen tensor (HPT) of the form  \emph{Spoligotypes} $\times$ \emph{Countries} $\times$ \emph{Time}. The final HPT is denoted as $\mathbf{\underline{X}} \in \mathbb{R}^{(I=311) \times (J=104) \times (K=7)}$. The host-pathogen tensor (HPT) is shown in Figure \ref{HPT}.

\begin{figure}[!]
\centering
\includegraphics[width=3.0in]{Figures/HPT.pdf}
\caption{Host-pathogen tensor (HPT). The first mode represents spoligotypes, the second mode represents countries, and the third mode represents time. This HPT of the form \emph{Spoligotypes} $\times$ \emph{Countries} $\times$ \emph{Time} is denoted as $\mathbf{\underline{X}} \in \mathbb{R}^{(I=311) \times (J=104) \times (K=7)}$.}
\label{HPT}
\end{figure}

\subsection{Distance matrices}
In the host-pathogen tensor, the first mode represents pathogen attributes, in this case spoligotypes. Genetic proximity of spoligotypes can be found using genetic distance measures. Hosts with close spoligotypes are more likely to be involved in the same mutation event. Similarly, the second mode represents host attributes, in this case country of birth. Proximity of countries can be found based on neighbourhood. Patients from close countries based on the proximity values are more likely to be involved in the same transmission event.

\subsubsection{Genetic proximity matrix}
Given 311 distinct spoligotypes, we define a genetic proximity measure between them. Mutation of spoligotypes is based on the Contiguous Deletion Assumption (CDA), which states that one or more contiguous spacers can be deleted in a mutation event, but not gained. Let $s_i$ represent spoligotype $i$, and let $s_i \rightarrow s_j$ represent the mutation of spoligotype $s_i$ into spoligotype $s_j$. Then, we define the CDA matrix, which summarizes contiguous deletion assumption, as follows:

\begin{equation*}\label{CDA}
%\boxed{
\text{CDA}(s_i, s_j) = \begin{cases}
\text{true} \text{,} & \quad \text{if } s_i \rightarrow s_j \text{ or } s_j \rightarrow s_i\\
\text{false} \text{,} & \quad \text{otherwise.} \end{cases}
%}
\end{equation*}

Let $H(s_i, s_j)$ be the Hamming distance between spoligotypes $s_i$ and $s_j$, as defined in \cite{InsightsIntoBIBM2011}:

\begin{equation*}
H\left(s_i, s_j\right) = \displaystyle\sum_{r=1}^{43} \mid s_{ir} - s_{jr}\mid
\end{equation*}

\noindent where $s_{ir}$ represents the value of $r-th$ spacer of spoligotype $s_i$. Then, we define the genetic proximity matrix $P_G$ as follows:

\begin{equation*}\label{P_G}
%\boxed{
P_G(s_i, s_j) = \begin{cases}
\displaystyle \frac{1}{1+H(s_i,s_j)} \text{,} & \quad \text{if } i\neq j \text{, } \text{CDA}(s_i, s_j) \text{, } H(s_i, s_j) \leq 10 \\
1 \text{,} & \quad \text{if } i=j \\
0 \text{,} & \quad \text{otherwise.} \end{cases}
%}
\end{equation*}

Genetic proximity matrix $P_G$ has values inversely proportional to the Hamming distance between two spoligotypes, as long as the Hamming distance between them is at most 10. For spoligotype pairs with $H(s_i, s_j) > 10$, the genetic proximity is set to zero. As a result, genetic proximity matrix reflects the likelihood of two different pathogens being involved in the same mutation event.

\subsubsection{Spatial proximity matrix}
Given 104 countries, we first define the Country Neighbourhood Matrix (CNM). Given two countries $C_i$ and $C_j$, the CNM is defined as follows:

%\begin{figure}[!]
\begin{equation*}\label{CNM}
%\boxed{
\text{CNM}(C_i, C_j) = \begin{cases}
\displaystyle 1 \text{,} & \quad \text{if } C_i \text{ and } C_j \text{ are neighbours} \\
0 \text{,} & \quad \text{otherwise.} \end{cases}
%}
\end{equation*}
%\end{figure}

Let $L(C_i, C_j)$ be the length of shortest path from $C_i$ to $C_j$ based on Dijkstra's shortest path algorithm on CNM. Then, we define the spatial proximity matrix $P_S$ as follows:

\begin{equation*}\label{P_S}
%\boxed{
P_S(C_i, C_j) = \begin{cases}
\displaystyle \frac{1}{1+L(C_i,C_j)} \text{,} & \quad \text{if } i\neq j \text{, } L(C_i, C_j) \leq 3 \\
1 \text{,} & \quad \text{if } i=j \\
0 \text{,} & \quad \text{otherwise.} \end{cases}
%}
\end{equation*}

Spatial proximity matrix $P_S$ has values inversely proportional to the length of shortest path between two countries, as long as the shortest-path length is at most 3. For country pairs with shortest-path length $L(C_i, C_j) > 3$, the proximity between two countries is set to zero. As a result, spatial proximity matrix reflects the likelihood of patients from two countries being involved in the same transmission event.

\begin{figure}[!]
\centering
\includegraphics[width=6.0in]{Figures/UBF.pdf}
\caption{Unified Biclustering Framework (UBF). In the first step, the data is generated as a matrix, a tensor, a coupled matrix-matrix, or a coupled matrix-tensor. In the second step, the data in various forms are factorized. In the third step, feature pattern similarity matrix is generated using the factor matrices of the decomposition. In the fourth step, we bicluster data points using \texttt{DensityInvariantBiclustering} algorithm. In the final step, we find the most stable biclusters using average best-match score.}
\label{UBF}
\end{figure}

\subsection{UBF: Unified Biclustering Framework}
In order to analyze host-pathogen associations using various forms of the raw dataset, we propose the Unified Biclustering Framework (UBF). Based on this framework, we generate the data in the first step, which can be a matrix, a tensor, a coupled matrix-matrix, or a coupled matrix-tensor. In the second step, we decompose datasets of various forms according to their type. In the third step, we generate the feature pattern similarity matrix. Fourth, we run the Density Invariant Biclustering (DIB) algorithm on the feature pattern similarity matrix. Finally, we find statistically significant biclusters and evaluate their biological relevance. Figure \ref{UBF} shows the steps of UBF. Next, we give the details of each step.

%The next two steps are borrowed from our earlier work, the Tensor Clustering Framework (TCF) \cite{TCF_MBT_BMC}. In the third step, we cluster the data points in the mode of interest using \texttt{kmeans\_mtimes\_seeded}, and validate clusters using Best-Match stability and F-measure.

\subsubsection{Data generation}
The host-pathogen tensor $\mathbf{\underline{X}} \in \mathbb{R}^{I \times J \times K}$ can be coupled with genetic proximity matrix $\mathbf{Y} \in \mathbb{R}^{I \times M}$ and spatial proximity matrix $\mathbf{Z} \in \mathbb{R}^{J \times N}$. This flexibility leads to different data configurations which allows simultaneous factorization of different data blocks. Possible data configurations are shown in Table \ref{DataConfigurations}. In data configuration 1, the host-pathogen tensor $\mathbf{\underline{X}}$ is summed and contracted along the time mode, and $\mathbf{\hat{X}}\in \mathbb{R}^{I \times J}$ is obtained and used without factorization. In data configuration 2, the original host-pathogen tensor $\mathbf{\underline{X}}$ is used. In configuration 3, genetic proximity matrix $\mathbf{Y}$ is coupled with the host-pathogen tensor $\mathbf{\underline{X}}$ in the first mode, incorporating the genetic distance to domain knowledge. In configuration 4, spatial proximity matrix $\mathbf{Z}$ is coupled with the host-pathogen tensor $\mathbf{\underline{X}}$ in the second mode, incorporating the spatial distance to domain knowledge. In configuration 5, genetic proximity matrix $\mathbf{Y}$ and spatial proximity matrix $\mathbf{Z}$ are coupled with the host-pathogen tensor $\mathbf{\underline{X}}$ in the first and second mode respectively, incorporating both the genetic distance and spatial distance to domain knowledge. In configuration 6, the host-pathogen tensor $\mathbf{\underline{X}}$ is contracted and summed along the time mode, keeping the genetic proximity matrix $\mathbf{Y}$ and spatial proximity matrix $\mathbf{Z}$ coupled with the contracted host-pathogen tensor, which is now the matrix $\mathbf{\hat{X}}$. We use all six configurations to find biclusters to associate spoligotypes and country of birth of tuberculosis patients, and test the effect of distance measures and time on detected groups.

%\begin{figure}[!t]
%\centering
%\includegraphics[width=3.5in]{Figures/DataConfigurationsWithBorder.pdf}
%\caption{Data configurations}
%\label{DataConfigurations}
%\end{figure}

% \raisebox{2ex - \height}

\begin{table}[!]
\centering
\begin{tabular}{|c|c|l|c|}
\cline{1-4}
Number & Data configuration & Extra information & Method in UBF \\ \hline \cline{1-4}
1 & \includegraphics[width=1.0in]{Figures/DC1_MBF.pdf} & --- & MBF \\ \cline{1-4}
2 & \includegraphics[width=1.0in]{Figures/DC2_TBF.pdf} & Time & TBF \\ \cline{1-4}
3 & \includegraphics[width=1.0in]{Figures/DC3_CMTBFg.pdf} & Time + genetic distance  & CMTBF$_{g}$ \\ \cline{1-4}
4 & \includegraphics[width=1.0in]{Figures/DC4_CMTBFs.pdf} & Time + spatial distance & CMTBF$_{s}$ \\ \cline{1-4}
5 & \includegraphics[width=1.0in]{Figures/DC5_CMTBFgs.pdf} & Time + genetic distance + spatial distance & CMTBF$_{gs}$ \\ \cline{1-4}
6 & \includegraphics[width=1.0in]{Figures/DC6_CMMBF.pdf} & Genetic distance + spatial distance & CMMBF \\ \cline{1-4}
\end{tabular}
\caption{Data configurations. The mode name S represents spoligotypes, C represents countries, and T represents time in years. The first configuration is a raw $Spoligotypes \times Countries$ matrix, and analyzed using Matrix Biclustering Framework (MBF) as part of UBF. The second data configuration includes time information as the third mode of the tensor, and analyzed using Tensor Biclustering Framework (TBF) as part of UBF. Third, fourth and fifth data configurations are the results of concatenating the genetic proximity matrix, spatial proximity matrix, and both respectively, to the host-pathogen tensor. They are analyzed using Coupled Matrix-Tensor Biclustering Framework (CMTBF) as part of UBF. Finally, in data configuration 6, we exclude time information and analyze the resulting data using coupled matrix-matrix biclustering framework (CMMBF) as part of UBF.}
\label{DataConfigurations}
\end{table}


\subsubsection{Data factorization}
In the second step of UBF, we factorize the dataset according to its form. If the data is a matrix, we use it as is. If it is a tensor, we use tensor decomposition methods, PARAFAC and Tucker3, and find the factor matrices for each mode. When the dataset is a coupled matrix-matrix or matrix-tensor, then we need to simultaneously factorize multiple matrices and/or tensors. We adopt the alternating least squares approach to solve coupled data factorizations. Next, we briefly outline the algorithms we use for coupled matrix-matrix factorization and coupled matrix-tensor factorization.

\paragraph{\textbf{Coupled matrix-matrix factorization (CMMF)}} Coupled matrices are simultaneously factorized using the CMMF\_ALS algorithm, which we outline next.
\subparagraph{\textbf{CMMF\_ALS}} The host-pathogen tensor contracted along the time mode becomes the matrix $\mathbf{\hat{X}} \in \mathbb{R}^{I \times J}$. Genetic proximity matrix $\mathbf{Y} \in \mathbb{R}^{I \times I}$ and spatial proximity matrix $\mathbf{Z} \in \mathbb{R}^{J \times J}$ are approximated as in the system of equations (\ref{ApproximationCMMF}).

\begin{align}\label{ApproximationCMMF}
\mathbf{\hat{X}} &\approx \mathbf{A} \mathbf{B}' & \nonumber\\
\mathbf{Y} &\approx \mathbf{A} \mathbf{V}' & \nonumber \\
\mathbf{Z} &\approx \mathbf{B} \mathbf{W}' &
\end{align}

\noindent We want to minimize the following loss function $L_1$, the sum of Frobenius norm of residuals for each data block:

\begin{align} \label{LossFunction}
L_1 = || \mathbf{\hat{X}} - \mathbf{A} \mathbf{B}'||_{F}^{2} +
     || \mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} +
     || \mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2} \nonumber \\
\end{align}

\noindent To minimize $L_1$, we first initialize the factor matrices $\mathbf{A},\mathbf{B},\mathbf{V},\mathbf{W}$ using truncated SVD, and then alternately minimize the loss function by fixing one of them at a time.

\begin{align}\label{CMMF}
\min_{A,B,V,W} & || \mathbf{\hat{X}} - \mathbf{A} \mathbf{B}'||_{F}^{2} + || \mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} + || \mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2} \nonumber \\
\min_{A,B,V,W} & \text{tr}\left( \left(\mathbf{\hat{X}} - \mathbf{A} \mathbf{B}'\right) \left( \mathbf{\hat{X}}' - \mathbf{B} \mathbf{A}'\right)\right)
               \text{tr}\left( \left(\mathbf{Y} - \mathbf{A} \mathbf{V}'\right) \left( \mathbf{Y}' - \mathbf{V} \mathbf{A}'\right)\right)
               \text{tr}\left( \left(\mathbf{Z} - \mathbf{B} \mathbf{W}'\right) \left( \mathbf{Z}' - \mathbf{W} \mathbf{B}'\right)\right) \nonumber \\
\min_{A,B,V,W} &
\text{tr}\left(\mathbf{\hat{X}} \mathbf{\hat{X}}'\right) -2\text{tr}\left(\mathbf{B} \mathbf{A}' \mathbf{\hat{X}}\right) +\text{tr}\left(\mathbf{A} \mathbf{B}' \mathbf{B} \mathbf{A}'\right) +
\text{tr}\left(\mathbf{Y} \mathbf{Y}'\right) -2\text{tr}\left(\mathbf{V} \mathbf{A}' \mathbf{Y}\right) +\text{tr}\left(\mathbf{A} \mathbf{V}' \mathbf{V} \mathbf{A}'\right) + \nonumber\\
&\text{tr}\left(\mathbf{Z} \mathbf{Z}'\right) -2\text{tr}\left(\mathbf{W} \mathbf{B}' \mathbf{Z}\right) +\text{tr}\left(\mathbf{B} \mathbf{W}' \mathbf{W} \mathbf{B}'\right) \nonumber \\
\min_{A,B,V,W} &
-2\text{tr}\left(\mathbf{B} \mathbf{A}' \mathbf{\hat{X}}\right) -2\text{tr}\left(\mathbf{V} \mathbf{A}' \mathbf{Y}\right) -2\text{tr}\left(\mathbf{W} \mathbf{B}' \mathbf{Z}\right) +
\text{tr}\left(\mathbf{A} \mathbf{B}' \mathbf{B} \mathbf{A}'\right) + \text{tr}\left(\mathbf{A} \mathbf{V}' \mathbf{V} \mathbf{A}'\right) +\text{tr}\left(\mathbf{B} \mathbf{W}' \mathbf{W} \mathbf{B}'\right)
\end{align}

\noindent The Lagrangian of the objective function (\ref{CMMF}) is:

\begin{equation}\label{LagrangianCMMF}
L = -2\text{tr}\left(\mathbf{B} \mathbf{A}' \mathbf{\hat{X}}\right) -2\text{tr}\left(\mathbf{V} \mathbf{A}' \mathbf{Y}\right) -2\text{tr}\left(\mathbf{W} \mathbf{B}' \mathbf{Z}\right) +
\text{tr}\left(\mathbf{A} \mathbf{B}' \mathbf{B} \mathbf{A}'\right) + \text{tr}\left(\mathbf{A} \mathbf{V}' \mathbf{V} \mathbf{A}'\right) +\text{tr}\left(\mathbf{B} \mathbf{W}' \mathbf{W} \mathbf{B}'\right)
\end{equation}

To minimize the loss function for $\mathbf{A},\mathbf{B},\mathbf{V},\mathbf{W}$ after fixing other factor matrices, we take the derivative of Lagrangian function $L$ in equation (\ref{LagrangianCMMF}), and set it to zero for each factor matrix, which gives the following update rules of matrices in CMMF\_ALS:

\begin{align*}
\mathbf{A} &= \left(\mathbf{X} \mathbf{B} + \mathbf{Y} \mathbf{V}\right) \backslash \left(\mathbf{B}' \mathbf{B} + \mathbf{V}' \mathbf{V} \right) \\
\mathbf{B} &= \left(\mathbf{X}' \mathbf{A} + \mathbf{Z} \mathbf{W}\right) \backslash \left(\mathbf{A}' \mathbf{A} + \mathbf{W}' \mathbf{W} \right) \\
\mathbf{V} &= \left(\mathbf{Y}' \mathbf{A}\right) \backslash \left(\mathbf{A}' \mathbf{A}\right) \\
\mathbf{W} &= \left(\mathbf{Z}' \mathbf{B}\right) \backslash \left(\mathbf{B}' \mathbf{B}\right)
\end{align*}

\noindent where $\backslash$ represents right matrix division. The complete CMMF\_ALS procedure is summarized in Algorithm \ref{CMMFALS}. In this algorithm, the function \texttt{svd\_mmf}($\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z}$) initializes the factor matrices $\mathbf{A},\mathbf{B},\mathbf{V},\mathbf{W}$ using truncated SVD with $min(J,M,N)$ components.

\begin{program}
%\begin{program}[H]
\caption{\small\texttt{CMMF\_ALS($\mathbf{X}\in \mathbb{R}^{I \times J}$,$\mathbf{Y}\in \mathbb{R}^{I \times M}$,$\mathbf{Z}\in \mathbb{R}^{J \times N}$)}}
\label{CMMFALS}
\begin{algorithmic}[1]
%\REQUIRE {\small\texttt{StrainDataset}} with spoligotypes and MIRU patterns.
%\ENSURE Spoligoforest $G = \left(V,E\right)$, where node set $V$ represents spoligotypes, and edge set $E$ represents spoligotype mutations.
\STATE [$\mathbf{A}$, $\mathbf{B}$, $\mathbf{V}$, $\mathbf{W}$] = \texttt{svd\_mmf}($\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z}$);
\STATE loss(current) = $|| \mathbf{X} - \mathbf{A} \mathbf{B}'||_{F}^{2} + || \mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} + || \mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2}$
\STATE $count=0$
\WHILE{(($count==0$) $\mid\mid$ ($0 < count \leq 10^3$ \&\& $\frac{|loss(current)-loss(prev)|}{loss(prev)} > 10^{-8}$))}
    \STATE $count++$
    \STATE // Solve for A
        \STATE $\mathbf{A} = \left(\mathbf{X} \mathbf{B} + \mathbf{Y} \mathbf{V}\right) \backslash \left(\mathbf{B}' \mathbf{B} + \mathbf{V}' \mathbf{V} \right)$
    \STATE // Solve for B
        \STATE $\mathbf{B} = \left(\mathbf{X}' \mathbf{A} + \mathbf{Z} \mathbf{W}\right) \backslash \left(\mathbf{A}' \mathbf{A} + \mathbf{W}' \mathbf{W} \right)$
    \STATE // Solve for V
        \STATE $\mathbf{V} = \left(\mathbf{Y}' \mathbf{A}\right) \backslash \left(\mathbf{A}' \mathbf{A}\right)$
    \STATE // Solve for W
        \STATE $\mathbf{W} = \left(\mathbf{Z}' \mathbf{B}\right) \backslash \left(\mathbf{B}' \mathbf{B}\right)$
    \STATE loss(current) = $|| \mathbf{X} - \mathbf{A} \mathbf{B}'||_{F}^{2} + || \mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} + || \mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2}$
\ENDWHILE
\end{algorithmic}
%\end{program}
\end{program}

\paragraph{\textbf{Coupled matrix-tensor factorization (CMTF)}} Coupled matrices and tensors can be simultaneously factorized. For this purpose, we used modifications of PARAFAC and Tucker3 methods. CMTF\_PARAFAC\_ALS decomposes the tensor using PARAFAC while factorizing the coupled matrices simultaneously. CMTF\_PARAFAC\_ALS algorithm and its variations exist in the literature. We build another algorithm, extension of Tucker3 to coupled matrix-tensor factorization. CMTF\_Tucker\_ALS algorithm decomposes the tensor using Tucker3, while simultaneously factorizing the coupled matrices. In the next section, we give the details of these algorithms.

\subparagraph{\textbf{CMTF\_PARAFAC\_ALS}}
Given the host-pathogen tensor $\mathbf{\underline{X}} \in \mathbb{R}^{I \times J \times K}$ coupled with genetic proximity matrix $\mathbf{Y} \in \mathbb{R}^{I \times I}$ and spatial proximity matrix $\mathbf{Z} \in \mathbb{R}^{J \times J}$, we approximate them as follows:

\begin{align}\label{ApproximationPARAFAC}
\mathbf{X_{(1)}} &\approx \mathbf{A} \left(\mathbf{C} \odot \mathbf{B}\right)'& \nonumber\\
\mathbf{Y} &\approx \mathbf{A} \mathbf{V}' & \nonumber \\
\mathbf{Z} &\approx \mathbf{B} \mathbf{W}' &
\end{align}

\noindent where $\odot$ denotes Khatri-Rao product. We want to minimize the following loss function which is the sum of squared Frobenius norm of each data block:

\begin{equation} \label{LossFunction}
L_2 = ||\mathbf{X_{(1)}} - \mathbf{A} \left(\mathbf{C} \odot \mathbf{B}\right)'||_{F}^{2} +
     ||\mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} +
     ||\mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2} \nonumber \\
%\min_{A,B,C,V,W} & L(\mathbf{\underline{X}}, \mathbf{Y}, \mathbf{Z}, \mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{V}, \mathbf{W})
\end{equation}

CMTF\_PARAFAC\_ALS is also known as CMTF\_ALS algorithm in the literature, which is detailed in earlier studies \cite{CMTF_OPT}. Therefore, we skip the details of the algorithm, and only focus on the update step for each factor matrix. Minimization for $\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{V},\mathbf{W}$ alternately returns the following updates at each step of CMTF\_PARAFAC\_ALS:

\noindent Update for $\mathbf{A}$:
\begin{align*}
\min_{A} & ||\mathbf{X_{(1)}} - \mathbf{A} \left(\mathbf{C} \odot \mathbf{B}\right)'||_{F}^{2} +
           ||\mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} \\
\min_{A} & ||\underbrace{[\mathbf{X_{(1)}} \,\, \mathbf{Y}]}_{T} - \mathbf{A} \underbrace{\left[\left(\mathbf{C} \odot \mathbf{B}\right)' \,\, \mathbf{V}'\right]}_{K}||_{F}^{2} \\
\Longrightarrow &\mathbf{A} = \left(\mathbf{T} \mathbf{K}'\right) / \left(\mathbf{K} \mathbf{K}'\right)
\end{align*}

\noindent Update for $\mathbf{B}$:
\begin{align*}
\min_{B} & ||\mathbf{X_{(2)}} - \mathbf{B} \left(\mathbf{C} \odot \mathbf{A}\right)'||_{F}^{2} +
           ||\mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2} \\
\min_{B} & ||\underbrace{[\mathbf{X_{(2)}} \,\, \mathbf{Z}]}_{T} - \mathbf{B} \underbrace{\left[\left(\mathbf{C} \odot \mathbf{A}\right)' \,\, \mathbf{W}'\right]}_{K}||_{F}^{2} \\
\Longrightarrow &\mathbf{B} = \left(\mathbf{T} \mathbf{K}'\right) / \left(\mathbf{K} \mathbf{K}'\right)
\end{align*}

\noindent Update for $\mathbf{C}$:
\begin{align*}
\min_{C} & ||\underbrace{\mathbf{X_{(3)}}}_{T} - \mathbf{C} \underbrace{\left(\mathbf{B} \odot \mathbf{A}\right)'}_{K}||_{F}^{2}\\
\Longrightarrow &\mathbf{C} = \left(\mathbf{T} \mathbf{K}'\right) / \left(\mathbf{K} \mathbf{K}'\right)
\end{align*}

\noindent Update for $\mathbf{V}$:
\begin{align*}
\min_{V} & ||\mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} \\
\Longrightarrow &\mathbf{V} = \left(\left(\mathbf{A}' \mathbf{A}\right) \backslash \left(\mathbf{A}' \mathbf{Y}\right)\right)'
\end{align*}

\noindent Update for $\mathbf{W}$:
\begin{align*}
\min_{W} & ||\mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2} \\
\Longrightarrow &\mathbf{W} = \left(\left(\mathbf{B}' \mathbf{B}\right) \backslash \left(\mathbf{B}' \mathbf{Z}\right)\right)'
\end{align*}

\subparagraph{\textbf{CMTF\_Tucker\_ALS}}
Next, we extend Tucker3 method to CMTF\_Tucker\_ALS for coupled matrix-tensor decomposition. This algorithm comes with the flexibility of factorizing the tensor using different number of components for each mode, while simultaneously factorizing the coupled matrices. The host-pathogen tensor $\mathbf{\underline{X}} \in \mathbb{R}^{I \times J \times K}$, genetic proximity matrix $\mathbf{Y} \in \mathbb{R}^{I \times I}$ and spatial proximity matrix $\mathbf{Z} \in \mathbb{R}^{J \times J}$ are approximated as in the system of equations (\ref{ApproximationTucker}).

\begin{align}\label{ApproximationTucker}
\mathbf{X_{(1)}} &\approx \mathbf{A} \mathbf{G_{(1)}} \left(\mathbf{C}' \otimes \mathbf{B}'\right)& \nonumber\\
\mathbf{Y} &\approx \mathbf{A} \mathbf{V}' & \nonumber \\
\mathbf{Z} &\approx \mathbf{B} \mathbf{W}' &
\end{align}

\noindent where $\otimes$ denotes Kronecker product. Note that in the Tucker3 model, the factor matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$ are orthogonal. Then, tensor $\mathbf{\underline{X}} \in \mathbb{R}^{I \times J \times K}$ can be decomposed using a $(P,Q,R)$-component Tucker3 model, while simultaneously factorizing $\mathbf{Y} \in \mathbb{R}^{I \times I}$ and $\mathbf{Z} \in \mathbb{R}^{J \times J}$ with the factor matrices of the shared mode. We want to minimize the loss function $L_3$ in equation (\ref{LossFunctionCMTFTuckerALS}), which is the sum of squared Frobenius norm of each data block.

\begin{equation} \label{LossFunctionCMTFTuckerALS}
L_3 = ||\mathbf{X_{(1)}} - \mathbf{A} \mathbf{G_{(1)}} \left( \mathbf{C}' \otimes \mathbf{B}' \right)||_{F}^{2} +
     ||\mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} +
     ||\mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2} \\
%\min_{A,B,C,V,W} & L(\mathbf{\underline{X}}, \mathbf{Y}, \mathbf{Z}, \mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{V}, \mathbf{W})
\end{equation}

\noindent To minimize $L_3$, we first initialize the factor matrices $\mathbf{A},\mathbf{B},\mathbf{V},\mathbf{W}$ using truncated SVD, and then alternately minimize the loss function for one of the variables at a time, while fixing the other variables. The following steps in equation (\ref{CMTFtuckerALSsteps}) reformulate the minimization of the loss function.

%\begin{figure*}
\begin{align}\label{CMTFtuckerALSsteps}
\min_{A} & ||\mathbf{X_{(1)}} - \mathbf{A} \mathbf{G_{(1)}} \left(\mathbf{C}' \otimes \mathbf{B}'\right)||_{F}^{2} +
           ||\mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} \nonumber \\
\min_{A} & ||\left[\mathbf{X_{(1)}} \,\, \mathbf{Y}\right] - \left[\mathbf{A} \mathbf{G_{(1)}} \left(\mathbf{C}' \otimes \mathbf{B}'\right) \,\, \mathbf{A}\mathbf{V}'\right]||_{F}^{2} \nonumber \\
\min_{A} & ||\left[\mathbf{X_{(1)}} \,\, \mathbf{Y}\right] - \left[\mathbf{A} \mathbf{A}' \underbrace{\mathbf{X_{(1)}} \left(\mathbf{C}\mathbf{C}' \otimes \mathbf{B}\mathbf{B}'\right)}_{\mathbf{M}} \,\, \mathbf{A}\mathbf{V}'\right]||_{F}^{2} \nonumber \\
\min_{A} & ||\left[\mathbf{X_{(1)}} \,\, \mathbf{Y}\right] - \left[\mathbf{A} \mathbf{A}' \mathbf{M} \,\, \mathbf{A}\mathbf{V}'\right]||_{F}^{2} \nonumber \\
\min_{A} & +\text{tr}\left(\left(\left[\mathbf{X_{(1)}} \,\, \mathbf{Y}\right] - \left[\mathbf{A} \mathbf{A}' \mathbf{M} \,\, \mathbf{A}\mathbf{V}'\right]\right) \left(\left[\mathbf{X_{(1)}} \,\, \mathbf{Y}\right]' - \left[\mathbf{A} \mathbf{A}' \mathbf{M} \,\, \mathbf{A}\mathbf{V}'\right]'\right) \right) \nonumber \\
\min_{A} & +\text{tr}\left( \left[\mathbf{X_{(1)}} \,\, \mathbf{Y}\right] \left[\mathbf{X_{(1)}} \,\, \mathbf{Y}\right]'\right)
           -2\text{tr}\left( \left[\mathbf{X_{(1)}} \,\, \mathbf{Y}\right] \left[\mathbf{M}' \mathbf{A} \mathbf{A}' \,\,\text{;} \,\, \mathbf{V}\mathbf{A}'\right]\right)
           +\text{tr}\left( \left[\mathbf{A} \mathbf{A}' \mathbf{M} \,\, \mathbf{A}\mathbf{V}'\right] \left[\mathbf{M} \mathbf{A} \mathbf{A}' \,\,\text{;} \,\, \mathbf{V}\mathbf{A}'\right]\right) \nonumber \\
\min_{A} & -2\text{tr}\left( \left[\mathbf{X_{(1)}} \,\, \mathbf{Y}\right] \left[\mathbf{M}' \mathbf{A} \mathbf{A}' \,\,\text{;} \,\, \mathbf{V}\mathbf{A}'\right]\right)
           +\text{tr}\left( \mathbf{A} \mathbf{A}' \mathbf{M} \mathbf{M}' \mathbf{A} \mathbf{A}' + \mathbf{A} \mathbf{V}' \mathbf{V} \mathbf{A}'\right) \nonumber \\
\min_{A} & -2\text{tr}\left( \mathbf{X_{(1)}}\mathbf{M}' \mathbf{A} \mathbf{A}' + \mathbf{Y} \mathbf{V}\mathbf{A}'\right)
           +\text{tr}\left( \mathbf{A} \mathbf{A}' \mathbf{M} \mathbf{M}' \mathbf{A} \mathbf{A}' + \mathbf{A} \mathbf{V}' \mathbf{V} \mathbf{A}'\right) \nonumber \\
\min_{A} & -2\text{tr}\left( \mathbf{X_{(1)}}\mathbf{M}' \mathbf{A} \mathbf{A}'\right) -2\text{tr}\left(\mathbf{Y} \mathbf{V}\mathbf{A}'\right)
           +\text{tr}\left( \mathbf{A} \mathbf{A}' \mathbf{M} \mathbf{M}' \mathbf{A} \mathbf{A}'\right) + +\text{tr}\left(\mathbf{A} \mathbf{V}' \mathbf{V} \mathbf{A}'\right) \nonumber \\
\min_{A} & -2\text{tr}\left( \mathbf{M}\mathbf{M}' \mathbf{A} \mathbf{A}'\right) -2\text{tr}\left(\mathbf{Y} \mathbf{Y}' \mathbf{A} \mathbf{A}'\right)
           +\text{tr}\left( \mathbf{A}' \mathbf{M} \mathbf{M}' \mathbf{A} \right) + +\text{tr}\left(\mathbf{A} \mathbf{A}' \mathbf{Y} \mathbf{Y}' \mathbf{A} \mathbf{A}' \right) \nonumber \\
\min_{A} & -2\text{tr}\left( \mathbf{A}' \mathbf{M}\mathbf{M}' \mathbf{A} \right) -2\text{tr}\left(\mathbf{A}' \mathbf{Y} \mathbf{Y}' \mathbf{A} \right)
           +\text{tr}\left( \mathbf{A}' \mathbf{M} \mathbf{M}' \mathbf{A} \right) + +\text{tr}\left(\mathbf{A}' \mathbf{Y} \mathbf{Y}' \mathbf{A} \right) \nonumber \\
\min_{A} & -\text{tr}\left( \mathbf{A}' \mathbf{M}\mathbf{M}' \mathbf{A} \right) -\text{tr}\left(\mathbf{A}' \mathbf{Y} \mathbf{Y}' \mathbf{A} \right) \\
\text{s.t} & \,\,\, \mathbf{A}'\mathbf{A} = \mathbf{I} \nonumber
%\Longrightarrow &\mathbf{A} = \left(\mathbf{T} \mathbf{K}'\right) / \left(\mathbf{K} \mathbf{K}'\right)
\end{align}
%\end{figure*}

\noindent The Lagrangian of this function is:

\begin{equation*}
L = -\text{tr}\left( \mathbf{A}' \mathbf{M}\mathbf{M}' \mathbf{A} \right) -\text{tr}\left(\mathbf{A}' \mathbf{Y} \mathbf{Y}' \mathbf{A} \right)
    + \lambda \left(\text{tr}\left(\mathbf{A}' \mathbf{A} - \mathbf{I} \right) \right)\\
\end{equation*}

\noindent The derivative of $L$ with respect to $\mathbf{A}$ set to zero returns the following equation:

\begin{equation*}
\left(\mathbf{M} \mathbf{M}' + \mathbf{Y} \mathbf{Y}'\right) \mathbf{A} = \lambda \mathbf{A}
\end{equation*}

\noindent which means that $\mathbf{A}$ is composed of first $I$ eigenvectors of $\left(\mathbf{M} \mathbf{M}' + \mathbf{Y} \mathbf{Y}'\right)$. We denote it as follows:

\begin{align*}
\mathbf{M} &= \mathbf{X_{(1)}} \left(\mathbf{C}\mathbf{C}' \otimes \mathbf{B}\mathbf{B}'\right) \\
\mathbf{A} &= \text{EVD}\left(\mathbf{M} \mathbf{M}' + \mathbf{Y} \mathbf{Y}', I \right)
\end{align*}

\noindent Similarly, for the other coupled mode, $\mathbf{B}$ is updated as follows:

\begin{align*}
\mathbf{M_2} &= \mathbf{X_{(2)}} \left(\mathbf{C}\mathbf{C}' \otimes \mathbf{A}\mathbf{A}'\right) \\
\mathbf{B} &= \text{EVD}\left(\mathbf{M_2} \mathbf{M_2}' + \mathbf{Z} \mathbf{Z}', J \right)
\end{align*}

\noindent And for the uncoupled mode, $\mathbf{C}$ is updated using SVD as in Tucker3-ALS algorithm:

\begin{align*}
\mathbf{M_3} &= \mathbf{X_{(3)}} \left(\mathbf{B}\mathbf{B}' \otimes \mathbf{A}\mathbf{A}'\right) \\
\mathbf{C} &= \text{SVD}\left(\mathbf{M_3}, K \right)
\end{align*}

The complete CMTF\_Tucker\_ALS procedure using these update rules is summarized in Algorithm \ref{CMTFtuckerALS}. Note that the function call \texttt{hosvd\_Tucker}($\mathbf{\underline{X}}$, [$P,Q,R$]) at the beginning of the algorithm initializes factor matrices via truncated SVD using $P,Q,R$ components respectively for each mode. The function \texttt{unfoldall($\mathbf{\underline{X}}$)} matricizes the tensor along each mode.

\begin{program}
%\begin{program}[H]
\caption{\small\texttt{CMTF\_Tucker\_ALS($\mathbf{\underline{X}}\in \mathbb{R}^{I \times J \times K}$,$\mathbf{Y}\in \mathbb{R}^{I \times M}$,$\mathbf{Z}\in \mathbb{R}^{J \times N}$, [$P,Q,R$])}}
\label{CMTFtuckerALS}
\begin{algorithmic}[1]
\STATE [$\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, $\mathbf{\underline{G}}$] = \texttt{hosvd\_Tucker}($\mathbf{\underline{X}}$, [$P,Q,R$]);
\STATE $\mathbf{V} = ((\mathbf{A}'\mathbf{A}) \backslash (\mathbf{A}'\mathbf{Y}))'$
\STATE $\mathbf{W} = ((\mathbf{B}'\mathbf{B}) \backslash (\mathbf{B}'\mathbf{Z}))'$
\STATE $[\mathbf{X_{(1)}},\mathbf{X_{(2)}},\mathbf{X_{(3)}}]$ = \texttt{unfoldall($\mathbf{\underline{X}}$)}
\STATE $[\mathbf{G_{(1)}},\mathbf{G_{(2)}},\mathbf{G_{(3)}}]$ = \texttt{unfoldall($\mathbf{\underline{G}}$)}
\STATE loss(current) = $||\mathbf{X_{(1)}} - \mathbf{A} \mathbf{G_{(1)}} \left(\mathbf{C}' \otimes \mathbf{B}' \right)||_{F}^{2} +
       ||\mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} + ||\mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2}$
\STATE $count=0$
\WHILE{(($count==0$) $\mid\mid$ ($0 < count \leq 10^3$  \&\& $\frac{|loss(current)-loss(prev)|}{loss(prev)} > 10^{-8}$))}
    \STATE $count++$
    \STATE // Solve for A
        \STATE $\mathbf{M_1} = \mathbf{X_{(1)}} \left(\mathbf{C}\mathbf{C}' \otimes \mathbf{B}\mathbf{B}'\right)$
        \STATE $\mathbf{A} = \text{EVD}\left(\mathbf{M_1} \mathbf{M_1}' + \mathbf{Y} \mathbf{Y}', I \right)$
    \STATE // Solve for B
        \STATE $\mathbf{M_2} = \mathbf{X_{(2)}} \left(\mathbf{C}\mathbf{C}' \otimes \mathbf{A}\mathbf{A}'\right)$
        \STATE $\mathbf{B} = \text{EVD}\left(\mathbf{M_2} \mathbf{M_2}' + \mathbf{Z} \mathbf{Z}', J \right)$
    \STATE // Solve for C
        \STATE $\mathbf{M_3} = \mathbf{X_{(3)}} \left(\mathbf{B}\mathbf{B}' \otimes \mathbf{A}\mathbf{A}'\right)$
        \STATE $\mathbf{C} = \text{SVD}\left(\mathbf{M_3}, K \right)$
    \STATE // Solve for V
        \STATE $\mathbf{V} = ((\mathbf{A}'\mathbf{A}) \backslash (\mathbf{A}'\mathbf{Y}))'$
    \STATE // Solve for W
        \STATE $\mathbf{W} = ((\mathbf{B}'\mathbf{B}) \backslash (\mathbf{B}'\mathbf{Z}))'$
    \STATE loss(current) = $||\mathbf{X_{(1)}} - \mathbf{A} \mathbf{G_{(1)}} \left(\mathbf{C}' \otimes \mathbf{B}' \right)||_{F}^{2} +
           ||\mathbf{Y} - \mathbf{A} \mathbf{V}'||_{F}^{2} + ||\mathbf{Z} - \mathbf{B} \mathbf{W}'||_{F}^{2}$
\ENDWHILE
\end{algorithmic}
%\end{program}
\end{program}

\subsubsection{Feature pattern similarity matrix generation}
We calculate the similarity of feature patterns of a spoligotype $s$ and country $c$ by calculating cosine similarity between feature pattern vectors of them. This is calculated in different ways for different forms of input data. If the input data is a matrix, then the matrix is used as the feature pattern similarity matrix (FPM). If the data is in tensor form, then feature pattern matrix is calculated for PARAFAC as follows. Assume that $R$-component PARAFAC model on the data matrix returns factor matrix $\mathbf{A}\in \mathbb{R}^{I \times R}$ for the first mode and factor matrix $\mathbf{B}\in \mathbb{R}^{J \times R}$ for the second mode. Then, we first normalize the rows of $\mathbf{A}$ and $\mathbf{B}$, and calculate the feature pattern matrix $FPM$ as follows:

\begin{equation}\label{PARAFAC_FPM}
\text{FPM}_{ij} = \begin{cases}
\displaystyle\frac{\mathbf{A}_{i.}\,\, \mathbf{B}_{j.}'}{||\mathbf{A}_{i.}||\,\, ||\mathbf{B}_{j.}||} \text{,} & \quad \text{if } N(i,j) > 0 \\
0 \text{,} & \quad \text{otherwise.}
\end{cases}
\end{equation}

where $N(i,j)$ represents the number of patients from country $i$ infected with strain $j$. This calculation is equivalent to cosine similarity of $i$-th feature vector of $\mathbf{A}$ and $j$-th feature vector of $\mathbf{B}$, only if there is at least one patient from country $i$ infected with strain $j$. Calculation of feature pattern matrix after applying Tucker3 model is slightly different. Assume that $(P,Q,R)$-component Tucker3 model on the data matrix returns factor matrix $\mathbf{A}\in \mathbb{R}^{I \times P}$ for the first mode, factor matrix $\mathbf{B}\in \mathbb{R}^{J \times Q}$ for the second mode, and the core tensor $\mathbf{\underline{G}}\in \mathbb{R}^{P \times Q \times R}$. First, we contract and sum the core tensor $\mathbf{\underline{G}}$ along the third mode and obtain $\mathbf{\hat{G}}$ matrix to calculate the level of interaction between the factors of $\mathbf{A}$ and $\mathbf{B}$. We normalize the rows of $\mathbf{A} \mathbf{\hat{G}}$ and $\mathbf{B}$. Finally, we calculate the feature pattern matrix as in equation (\ref{Tucker_FPM}).

\begin{align}\label{Tucker_FPM}
\mathbf{\hat{G}}_{pq} &= \displaystyle \sum_{r=1}^{R}\mathbf{\underline{G}}_{pqr} \nonumber\\
\text{FPM}_{ij} &=
\begin{cases}
\displaystyle\frac{\mathbf{A}_{i.}\,\mathbf{\hat{G}}}{||\mathbf{A}_{i.}\,\mathbf{\hat{G}}||} \,\, \displaystyle\frac{\mathbf{B}_{j.}'}{||\mathbf{B}_{j.}||} \text{,} & \quad \text{if } N(i,j) > 0 \\
0 \text{,} & \quad \text{otherwise.}
\end{cases}
\end{align}

For coupled factorizations, we use the same equations. After coupled matrix-matrix decomposition, we use Equation (\ref{PARAFAC_FPM}) to find the feature pattern similarity matrix. For coupled matrix-tensor factorization, if CMTF\_PARAFAC\_ALS is used for factorization, then FPM is calculated using Equation (\ref{PARAFAC_FPM}). If CMTF\_Tucker\_ALS is used for factorization, then Equation (\ref{Tucker_FPM}) is used to calculate FPM.

\subsubsection{Density-Invariant Biclustering}
In this section, we introduce a novel biclustering algorithm based on an existing algorithm and several graph attributes. First, we discretize the input matrix and use it as input to BiMax algorithm to find inclusion-maximal biclusters \cite{BiMax}. Then, we use these biclusters as seed, and find density and variance of these biclusters, which are bicliques. Finally, we find density-invariant biclusters among candidate inclusion-maximal biclusters.

%\paragraph{Preprocessing}
%Given the feature pattern matrix $\mathbf{X} \in \mathbb{R}^{I \times J}$, we autoscale it to have zero mean, and variance of 1. This results in the following matrix $\mathbf{\bar{X}}$:
%
%\begin{equation*}
%\mathbf{\bar{X}}_{gc} = \frac{\mathbf{X}_{gc}-\mu_g}{\sigma_g}
%\end{equation*}
%
%\noindent where $\mu_g$ is the mean of $g$-th row of $\mathbf{X}$ associated with gene $g$, and $\sigma_g$ is the variance of $g$-th row of $\mathbf{X}$ associated with gene $g$.
%
Given the feature pattern similarity matrix $\mathbf{X} \in \mathbb{R}^{I \times J}$, we use density-invariant biclustering to find coherent biclusters. Let $G=(U, V, E)$ represent a bipartite graph, where $U$ represents the set of genes, or rows in $\mathbf{X}$, $V$ represents the set of conditions, or columns in $\mathbf{X}$, and $E$ represents the weight of the edges connecting vertex set $U$ and vertex set $V$. The weights $E$ are equivalent to values of matrix $\mathbf{X}$. We want to find biclusters of the following form:

\begin{equation}
B_i = (U_i, V_i, E_i)
\end{equation}

\noindent where $\displaystyle B=\bigcup_{i=1}^{n}B_i$ is a biclustering of rows and columns of $\mathbf{X}$. Each bicluster associates a set of rows, in this case spoligotypes, to a set of columns, in this case countries. Notice that each bicluster maps to a submatrix of the original data matrix.

Density-Invariant biclustering algorithm first discretizes edge weights using a weight threshold $th$, and converts the input matrix into a binary matrix $\mathbf{D}$. Then we use the binary inclusion-maximal biclustering algorithm (BiMax) by Prelic et al. on this binary matrix and find a set of candidate biclusters \cite{BiMax}. These biclusters are inclusion-maximal, because the submatrices corresponding to these biclusters are all 1's, and there is no other bicluster which is a superset of it. Output of BiMax algorithm after discretization returns a good starting point for density-invariant biclustering algorithm. Next, we focus on these candidate biclusters. For this purpose, we define density and variance of a graph.

\begin{mydef}
\textbf{Density of a graph:} Density of a graph is the average weight of its edges. Given a graph $G=(V,E)$ where w(e) represents the weight of edge $e \in E$, the density of graph $G$ is calculated as follows:
\begin{equation*}
d(G) = \frac{\displaystyle \sum_{e\in E}w(e)}{\binom{|V|}{2}}
\end{equation*}
\end{mydef}

\begin{mydef}
\textbf{Variance of a graph:} Variance of a graph is the standard deviation of its edge weights. Given a graph $G=(V,E)$ where w(e) represents the weight of edge $e \in E$, the variance of graph $G$ is calculated as follows:
\begin{equation*}
v(G) = \sqrt{\frac{1}{|E|-1} \displaystyle \sum_{e\in E}{\left( w\left(e\right) - \bar{w}\right)}^2}
\end{equation*}
\end{mydef}

Using the density and variance of graphs, we can define a new set of graphs which are bounded by their edge weights. We define $\alpha$-dense $\beta$-variant biclusters, or density-invariant biclusters, which are graphs of the form $B=(U,V,E)$ with density $d(B)\geq\alpha$ and variance $v(B)\leq\beta$, and similarly for all one-vertex-induced subgraphs of the $B=(U,V,E)$.

\begin{mydef}
\textbf{Density-invariant bicluster:} Let $B=(U,V,E)$ be a bicluster, where edges in $E$ connect vertices in $U$ to vertices in $V$. Bicluster $B$ is an $\alpha$-dense bicluster if $d(B) \geq \alpha$, and it is a $\beta$-variant bicluster if $v(B) \leq \beta$. Define $B'$ as an induced subgraph of $B$ after removing one vertex from $B$, either from vertex set $U$ or vertex set $V$. Bicluster $B=(U,V,E)$ is an ($\alpha,\beta$)-density-invariant bicluster, or density-invariant bicluster, if $B$ and all its one-vertex-induced subgraphs are $\alpha$-dense $\beta$-variant. In short, bicluster $B$ is a density-invariant bicluster if the following conditions hold:
\begin{enumerate}
\item $d(B) \geq \alpha$, $v(B) \leq \beta$
\item $d(B') \geq \alpha$, $v(B') \leq \beta$ $ \quad \forall B' = B\setminus\{m\}$ \text{where} $m \in U \cup V$, $|B'|>0$
\end{enumerate}
\end{mydef}

\noindent Notice that a density-invariant bicluster forms a biclique with average weight bounded from below, and variance of weights bounded from above. All induced subgraphs obtained after removing one vertex from a density-invariant bicluster are still $\alpha$-dense and $\beta$-variant, but not necessarily density-invariant biclusters. At this point, we define strong antimonotonicity of a graph, which was introduced in Pao et al. \cite{DCB_PaoRecep}.

\begin{mydef}
\textbf{Strong antimonotonicity:} A graph attribute is strong antimonotone if for each graph $G=(V,E)$ with the attribute, every induced subgraph $G' = G-\{v\}$, has the attribute, where $v \in V$.
\end{mydef}

According to the definition of strong antimonotonicity, the attribute of being a density-invariant graph or bicluster is not strongly antimonotone. This is because the vertex-induced subgraphs of the original graph are $\alpha$-dense and $\beta$-variant, but their vertex-induced subgraphs need not be ($\alpha,\beta$)-density-invariant biclusters.

Finally, we iterate over candidate biclusters found as output from BiMax algorithm and find density-invariant biclusters among these candidate biclusters. This results in strongly connected and more homogeneous biclusters. Algorithm \ref{DensityInvariantBiclustering} summarizes \texttt{DensityInvariantBiclustering} procedure.

\begin{program}
%\begin{program}[H]
\caption{\small \texttt{Biclusters} = \texttt{DensityInvariantBiclustering($\mathbf{X}\in \mathbb{R}^{I \times J}$, $th$, $\alpha$, $\beta$)}}
\label{DensityInvariantBiclustering}
\begin{algorithmic}[1]
\REQUIRE {\small Data matrix $\mathbf{X}\in \mathbb{R}^{I \times J}$, discretization threshold $th$, density threshold $\alpha$, variance threshold $\beta$.}
\ENSURE Density-invariant biclusters \texttt{Biclusters}.
\STATE $\mathbf{D} = \texttt{discretize}(\mathbf{X}, th)$
\STATE \texttt{CandidateBiclusters} = \texttt{BiMax}($\mathbf{D}$)
\STATE \texttt{Biclusters} = $\emptyset$
\FOR{i=1:1:length(\texttt{CandidateBiclusters})}
    \STATE $B(U,V,E)$ = \texttt{CandidateBiclusters}(i)
    \STATE $check1 = \left(d(B) \geq \alpha\right)$ \&\& $\left(v(B) \leq \beta\right)$
    \STATE $check2 = true$
    \STATE M = $U \cup V$
    \FOR{j=1:1:length(M)}
        \STATE $m = M(j)$
        \STATE $B' = B \setminus \{m\}$
        \IF{$\left(\left(B' \neq \emptyset\right) \,\,\, \&\& \,\,\, !\left( d(B') \geq \alpha \,\,\, \&\& \,\,\, v(B') \leq \beta \right)\right)$}
            \STATE $check2 = false$
            \STATE break
        \ENDIF
    \ENDFOR
    \IF{($check1$ \&\& $check2$)}
        \STATE \texttt{Biclusters} = \texttt{Biclusters} $\cup \{B\}$
    \ENDIF
\ENDFOR
\end{algorithmic}
%\end{program}
\end{program}

In \texttt{DensityInvariantBiclustering} algorithm, $\texttt{discretize}(\mathbf{X}, th)$ function discretizes the input data matrix as follows:

\begin{equation*}
\mathbf{D}_{ij} = \begin{cases}
1 \text{,} & \quad \text{if } \mathbf{X}_{ij} \geq th \\
0 \text{,} & \quad \text{otherwise.} \end{cases}
\end{equation*}

\texttt{BiMax} algorithm is run on this binary matrix $\mathbf{D}$, and inclusion-maximal biclusters are obtained. Then, among these candidate biclusters, density-invariant biclusters are found.

\subsubsection{Statistically significant bicluster selection}
In order to find statistically significant biclusters, we sample 90\% of the patients, and rerun the biclustering algorithm, and obtain 20 new biclusterings. Then, we calculate the stability of each density-invariant bicluster found in the previous step using average best-match score. First, we calculate the match score of two biclusters $B_1 = (G_1, C_1)$, $B_2 = (G_2, C_2)$, where $G_1$, $G_2$ represent gene sets and $C_1$, $C_2$ represent condition sets. Similar to Prelic et al. and Lie et al. \cite{BiMax, MatchScoreLiuWang}, the match score of biclusters $B_1 = (G_1, C_1)$, $B_2 = (G_2, C_2)$ is calculated as follows:

\begin{equation}
\texttt{match}(B_1, B_2) = \frac{|G_1 \cap G_2| + |C_1 \cap C_2|}{|G_1 \cup G_2| + |C_1 \cup C_2|}
\end{equation}

Let $M = \displaystyle\bigcup_{i=1}^{k}B^{*}_i$ be a biclustering of the subsample of the dataset. We compare a bicluster $B=(G,C)$ to all biclusters in $B^{*}_i \in M$, and assign the maximum match value as the best-match score:

\begin{equation}
\texttt{best\_match}(B, M = \displaystyle\bigcup_{i=1}^{k}B^{*}_i) = \max_{B^{*}_i \in M} \texttt{match}(B, B^{*}_i)
\end{equation}

Finally, we take the average best-match score of each bicluster $B$ by comparing them to each biclustering $M_i$, and obtain the average best-match score of bicluster $B$ as follows:

\begin{equation}
\texttt{average\_best\_match}(B, \displaystyle\bigcup_{i=1}^{n}M_i) = \frac{\displaystyle\sum_{i=1}^{n}\texttt{best\_match}(B, M_i)}{n}
\end{equation}

Finally, we pick the biclusters with $\geq 95\%$ average best-match score as statistically significant biclusters, and evaluate their biological relevance. If there are no significant biclusters, we report top 5 stable biclusters with their average best-match scores.

\section{Results}
In order to find existing or hidden host-pathogen associations, we biclustered spoligotypes and countries using six different data configurations shown in Figure \ref{DataConfigurations}. For each data configuration, we followed the steps of Unified Biclustering Framework (UBF), and found the most stable biclusters. Table \ref{BiclusteringResults} shows the parameters of \texttt{DensityInvariantBiclustering} ($th$, $\alpha$, $\beta$) and number of density-invariant biclusters for each data configuration. Note that PARAFAC and Tucker3 variants of TBF, CMTF\_PARAFAC\_ALS and CMTF\_Tucker\_ALS variants of CMTBF are evaluated separately. Next, we evaluate the statistical significance and biological relevance of biclusters for each data configuration, and find host-pathogen associations within the whole patient dataset and within each major lineage.

\begin{table}[!]
\centering
\begin{tabular}{|c|l|c|c|}
\cline{1-4}
\# Configuration & Method & DIB parameters ($th$, $\alpha$, $\beta$) & \# DIB \\ \hline \hline
                   1 & MBF & 0.80, 0.80, 0.15 & 8 \\ \cline{1-4}
{\multirow{2}{*}{2}} & TBF (PARAFAC) & 0.98, 0.89, 0.01 & 170 \\ \cline{2-4}
                     & TBF (Tucker3) & 0.60, 0.60, 0.40 & 5 \\ \cline{1-4}
{\multirow{2}{*}{3}} & CMTBF$_g$ (CMTF\_PARAFAC\_ALS) & 0.60, 0.60, 0.40 & 0 \\ \cline{2-4}
                     & CMTBF$_g$ (CMTF\_Tucker\_ALS) & 0.70, 0.70, 0.30 & 4 \\ \cline{1-4}
{\multirow{2}{*}{4}} & CMTBF$_s$ (CMTF\_PARAFAC\_ALS) & 0.80, 0.90, 0.10 & 21 \\ \cline{2-4}
                     & CMTBF$_s$ (CMTF\_Tucker\_ALS) & 0.80, 0.85, 0.15 & 6 \\ \cline{1-4}
{\multirow{2}{*}{5}} & CMTBF$_{gs}$ (CMTF\_PARAFAC\_ALS) & 0.98, 0.99, 0.01 & 0 \\ \cline{2-4}
                     & CMTBF$_{gs}$ (CMTF\_Tucker\_ALS) & 0.60, 0.70,0.30 & 5 \\ \cline{1-4}
                   6 & CMMBF & 0.60, 0.60, 0.40 & 17 \\ \cline{1-4}
\end{tabular}
\caption{Biclustering results for each data configuration, including \texttt{DensityInvariantBiclustering} parameters and number of density-invariant biclusters (DIB). For TBF, PARAFAC and Tucker3 model results are listed separately. For CMTBF, CMTF\_PARAFAC\_ALS and CMTF\_Tucker\_ALS results are listed separately. When there are no stable biclusters with average best-match score $\geq 95\%$, five most stable biclusters are picked as the stable biclusters.}
\label{BiclusteringResults}
\end{table}

\subsection{Biclusters using spoligotypes and country of birth}
We first contract and sum the host-pathogen tensor along the time mode and find biclusters based on the distribution of spoligotypes to countries of birth, as in data configuration 1 in Figure \ref{DataConfigurations}. In this setting, no distance measure or time is added to the domain knowledge. Table \ref{BiclusteringResults_1} shows the density-invariant biclusters. Bicluster B1 suggests that patients from Haiti are infected with ST1162 strain, a Beijing strain, and ST398, a LAM4 strain. Bicluster B12 is listed in the supplementary material due to its size: \url{http://tbinsight.cs.rpi.edu/UBFsupp.rar}. This bicluster contains 848 patients from United States who are infected with 63 different strains. One of these strains is the transmissive Beijing strain ST1 which initiated many outbreaks in United States \cite{W-Beijing-GlobalDissemination, W-Beijing-Review}. Bicluster B13 shows that patients from China are infected with 7 different strains. Bicluster B14 shows that ST1908 and ST58 are two Ecuadorian isolates belonging to Euro-American lineage. Bicluster B16 is a well-known association, and suggests that patients from Philippines are infected with an EAI2-Manila strain, ST897. Bicluster B18 suggests that Mexican patients, as neighbours of United States, are infected with UST251 and ST1154, two Euro-American strains. The five most stable biclusters are B16, B17, B18, B11, B14, and their average best-match scores are in the range [0.1667, 0.2]. One may argue that biclusters with few patients does not constitute a strong host-pathogen association. This suggests that TB detection rate should be increased to gather more patient data and make more accurate inferences on host-pathogen association.

\begin{table}[!]
\centering
\begin{tabular}{|c|c||c|c|c||c|c|}
\cline{1-7}
{\multirow{2}{*}{\textbf{Bicluster}}} & {\multirow{1}{*}{\textbf{Number}}} & \multicolumn{3}{|c||}{\textbf{Spoligotypes}} & \multicolumn{2}{|c|}{\textbf{Countries}} \\ \cline{3-7}
                                      & {\multirow{1}{*}{\textbf{of patients}}} & \textbf{SIT no} & \textbf{Major lineage} & \textbf{Sublineage} & \textbf{Name} & \textbf{TB continent} \\ \cline{1-7}
{\multirow{2}{*}{B11}} & {\multirow{2}{*}{5}} & ST1162 & East-Asian & Beijing & {\multirow{2}{*}{Haiti}} & {\multirow{2}{*}{Americas}} \\ %\cline{2-6}
                       &                       & ST398 & Euro-American & LAM4 & & \\ \cline{1-7}
{\multirow{7}{*}{B13}} & {\multirow{7}{*}{19}} & ST265 & East-Asian & Beijing & {\multirow{7}{*}{China}} & {\multirow{7}{*}{East Asia}} \\ %\cline{2-6}
                       &                       & ST422 & \emph{M. bovis} & BOV\_1 & & \\
                       &                       & ST89 & Indo-Oceanic & EAI5 & & \\
                       &                       & ST287 & Indo-Oceanic & EAI2-Manila & & \\
                       &                       & ST1268 & Euro-American & T5 & & \\
                       &                       & ST25 & East-African Indian & CAS1-Delhi & & \\
                       &                       & ST732 & Euro-American & T1 & & \\
                       \cline{1-7}
{\multirow{2}{*}{B14}} & {\multirow{2}{*}{6}} & ST1908 & Euro-American & H3 & {\multirow{2}{*}{Ecuador}} & {\multirow{2}{*}{Americas}} \\ %\cline{2-6}
                       &                      & ST58 & Euro-American & T5 & & \\ \cline{1-7}
{\multirow{3}{*}{B15}} & {\multirow{3}{*}{6}} & ST43 & Indo-Oceanic & EAI6-BGD1 & {\multirow{1}{*}{Dominican}} & {\multirow{2}{*}{Americas}} \\ %\cline{2-6}
                       &                      & ST848 & Euro-American & T2 & {\multirow{1}{*}{Republic}} & \\
                       &                      & ST511 & Euro-American & H3 & & \\
                       \cline{1-7}
{\multirow{1}{*}{B16}} & {\multirow{1}{*}{2}} & ST897 & Indo-Oceanic & EAI2-Manila & {\multirow{1}{*}{Philippines}} & {\multirow{1}{*}{Southeast Asia}} \\ %\cline{2-6}
                       \cline{1-7}
{\multirow{1}{*}{B17}} & {\multirow{1}{*}{2}} & ST447 & Euro-American & T1 & {\multirow{1}{*}{Bangladesh}} & {\multirow{1}{*}{Indian Subcontinent}} \\ %\cline{2-6}
                       \cline{1-7}
{\multirow{2}{*}{B18}} & {\multirow{2}{*}{4}} & UST251 & Euro-American & S & {\multirow{2}{*}{Mexico}} & {\multirow{2}{*}{Americas}} \\ %\cline{2-6}
                       &                      & ST1154 & Euro-American & LAM9 & & \\
                       \cline{1-7}
\end{tabular}
\caption{Biclustering results on data configuration 1 using UBF. Biclusters associate spoligotypes to country of birth of patients. For spoligotypes, SIT number, major lineage based on CBN, and sublineage based on KBBN are listed. For countries, the name and the TB continent are listed. Biclusters B16 represents the well-known association between patients from Philippines and EAI2-Manila strains.}
\label{BiclusteringResults_1}
\end{table}

\subsection{Incorporating time}
The original host-pathogen tensor has time as the third mode. Therefore, when we found biclusters using the host-pathogen tensor as in data configuration 2 of Figure \ref{DataConfigurations}, we account for distribution of spoligotypes to countries of birth through time, in this case years from 2001 to 2007. When we use PARAFAC to decompose the host-pathogen tensor, we found 170 density-invariant biclusters. Here, we focus on five most stable biclusters when PARAFAC model is used. Average best-match scores of these five biclusters range from 0.6915 to 0.7295. The full list of these biclusters can be found in the supplementary material. Bicluster B211 associates Vietnamese patients to 11 strains belonging to Euro-American, East Asian, Indo-Oceanic and East-African Indian lineages. Bicluster B212 suggests that patients from Peru are infected with 17 different strains belonging to Euro-American, Indo-Oceanic, and East Asian lineages. Bicluster B214 is shown in Table \ref{BiclusteringResults_2}. There are 111 patients in bicluster B214 from India, Peru and Vietnam, which are infected with 6 Euro-American strains and one East Asian strain. Notice that this East Asian strain is ST1, which is the transmissive Beijing strain. This suggests that some of the patients in this bicluster must be involved in the outbreaks in United States initiated by ST1 Beijing strains.

When Tucker3 model is used to decompose the host-pathogen tensor, we find 5 density-invariant biclusters. Their average-best match scores range from 0.04 to 0.18, which shows that biclusters found using Tucker3 model are less stable compared to the ones found using PARAFAC model. These five biclusters, bicluster B221 to B225, are listed in the supplementary material. Bicluster B221 suggest that US patients are infected with 31 different strains, and bicluster B222 suggests that Chinese patients are infected with 11 strains belonging to Euro-American, East-African Indian, and \emph{M. bovis} lineages. Note that no East Asian strain is associated with Chinese patients, which suggests that our biclustering analysis on the host-pathogen tensor is introducing noise when time is added into domain knowledge. Bicluster B225, also listed in Table \ref{BiclusteringResults_2}, has 4 patients and suggests that Mexican patients are infected with ST294, ST290 and  ST176 strains, all members of Euro-American lineage.

\begin{table}[!]
\centering
\begin{tabular}{|c|c||c|c|c||c|c|}
\cline{1-7}
{\multirow{2}{*}{\textbf{Bicluster}}} & {\multirow{1}{*}{\textbf{Number}}} & \multicolumn{3}{|c||}{\textbf{Spoligotypes}} & \multicolumn{2}{|c|}{\textbf{Countries}} \\ \cline{3-7}
                                      & {\multirow{1}{*}{\textbf{of patients}}} & \textbf{SIT no} & \textbf{Major lineage} & \textbf{Sublineage} & \textbf{Name} & \textbf{TB continent} \\ \cline{1-7}
{\multirow{2}{*}{B214}}& {\multirow{2}{*}{111}} & ST53 & Euro-American & T1 & {\multirow{1}{*}{India}} & {\multirow{1}{*}{Indian Subcontinent}} \\ %\cline{2-6}
                       &                        & ST17 & Euro-American & LAM2 & {\multirow{1}{*}{Peru}} & {\multirow{1}{*}{Americas}} \\
                       &                        & ST1 & East Asian & Beijing & {\multirow{1}{*}{Vietnam}} & {\multirow{1}{*}{Southeast Asia}} \\
                       &                        & ST197 & Euro-American & X3 & & \\
                       &                        & ST61 & Euro-American & LAM10-CAM & & \\
                       &                        & ST119 & Euro-American & X1 & & \\
                       &                        & ST42 & Euro-American & LAM9 & & \\
                       \cline{1-7}
{\multirow{3}{*}{B225}}& {\multirow{3}{*}{4}}   & ST294 & Euro-American & H3 & {\multirow{3}{*}{Mexico}} & {\multirow{3}{*}{Americas}} \\ %\cline{2-6}
                       &                        & ST290 & Euro-American & LAM9 & & \\
                       &                        & ST176 & Euro-American & LAM6 & & \\
                       \cline{1-7}
\end{tabular}
\caption{Biclustering results on data configuration 2 using PARAFAC and Tucker3 models via UBF on the host-pathogen tensor. Bicluster B214 associates patients from India, Peru and Vietnam to 6 Euro-American strains, and the transmissive East Asian Beijing strain ST1. Bicluster B224 groups Mexican patients infected with three different Euro-American strains.}
\label{BiclusteringResults_2}
\end{table}

\subsection{Incorporating time and distance measures}
Next, we concatenate distance matrices one at a time, and finally both of them, to the host-pathogen tensor. Concatenation of genetic proximity matrix results in data configuration 3, concatenation of spatial proximity matrix results in data configuration 4, and concatenation of both results in data configuration 5. We factorize these matrices using coupled matrix-tensor factorization via CMTF\_PARAFAC\_ALS and CMTF\_Tucker\_ALS, and report statistically significant and biologically relevant biclusters. The full list of biclusters can be found in the supplementary material.

If we use genetic distance matrix only, factorization via CMTF\_PARAFAC\_ALS results in no density-invariant biclusters. When the coupled matrix-tensor is decomposed via CMTF\_Tucker\_ALS, 4 stable clusters are found, the stability of which range from 0.08 to 0.19. Two of these biclusters, B321 and B323, are listed in Table \ref{BiclusteringResults_3_4_5}. Bicluster B321 groups 32 patients from Ecuador, infected with ST53, ST62, ST51, ST1908, which are all Euro-American strains. Notice that ST53 and ST51 are T1 strains, which is a class of not-well-defined Euro-American strains. Bicluster B323 contains 4 patients from Mexico, all infected with ST52, a Euro-American T2 strain.

\begin{table}[!]
\centering
\begin{tabular}{|c|c||c|c|c||c|c|}
\cline{1-7}
{\multirow{2}{*}{\textbf{Bicluster}}} & {\multirow{1}{*}{\textbf{Number}}} & \multicolumn{3}{|c||}{\textbf{Spoligotypes}} & \multicolumn{2}{|c|}{\textbf{Countries}} \\ \cline{3-7}
                                      & {\multirow{1}{*}{\textbf{of patients}}} & \textbf{SIT no} & \textbf{Major lineage} & \textbf{Sublineage} & \textbf{Name} & \textbf{TB continent} \\ \cline{1-7}
{\multirow{4}{*}{B321}}& {\multirow{4}{*}{32}}  & ST53 & Euro-American & T1 & {\multirow{4}{*}{Ecuador}} & {\multirow{4}{*}{Americas}} \\ %\cline{2-6}
                       &                        & ST62 & Euro-American & H1 & & \\
                       &                        & ST51 & Euro-American & T1 & & \\
                       &                        & ST1908 & Euro-American & H3 & & \\
                       \cline{1-7}
{\multirow{1}{*}{B323}}& {\multirow{1}{*}{4}}   & ST52 & Euro-American & T1 & {\multirow{1}{*}{Mexico}} & {\multirow{1}{*}{Americas}} \\ %\cline{2-6}
                       \cline{1-7}
{\multirow{4}{*}{B411}}& {\multirow{4}{*}{6}}   & {\multirow{4}{*}{ST908}} & {\multirow{4}{*}{Euro-American}} & {\multirow{4}{*}{LAM2}} & Dominican Rep. & Americas \\ %\cline{2-6}
                       &                        &                          &                                  &                         & Puerto Rico    & Americas \\ %\cline{2-6}
                       &                        &                          &                                  &                         & Trinidad and Tobago    & Americas \\ %\cline{2-6}
                       &                        &                          &                                  &                         & United States  & Americas \\ %\cline{2-6}
                       \cline{1-7}
{\multirow{4}{*}{B412}}& {\multirow{4}{*}{6}}   & {\multirow{4}{*}{ST904}} & {\multirow{4}{*}{Euro-American}} & {\multirow{4}{*}{T5}}   & Ecuador & Americas \\ %\cline{2-6}
                       &                        &                          &                                  &                         & Haiti    & Americas \\ %\cline{2-6}
                       &                        &                          &                                  &                         & Trinidad and Tobago    & Americas \\ %\cline{2-6}
                       &                        &                          &                                  &                         & United States  & Americas \\ %\cline{2-6}
                       \cline{1-7}
{\multirow{2}{*}{B414}}& {\multirow{2}{*}{6}}   & {\multirow{1}{*}{ST904}} & {\multirow{1}{*}{Euro-American}} & {\multirow{1}{*}{T5}}   & Trinidad and Tobago    & Americas \\ %\cline{2-6}
                       &                        & {\multirow{1}{*}{ST908}} & {\multirow{1}{*}{Euro-American}} & {\multirow{1}{*}{LAM2}} & United States          & Americas \\ %\cline{2-6}
                       \cline{1-7}
{\multirow{6}{*}{B421}}& {\multirow{6}{*}{32}}  & {\multirow{6}{*}{ST1}}   & {\multirow{6}{*}{East Asian}}    & {\multirow{6}{*}{Beijing}} & Taiwan      & East Asian \\ %\cline{2-6}
                       &                        &                          &                                  &                            & Barbados    & Americas \\ %\cline{2-6}
                       &                        &                          &                                  &                            & Dominica    & Americas \\ %\cline{2-6}
                       &                        &                          &                                  &                            & Malaysia    & Southeast Asia \\ %\cline{2-6}
                       &                        &                          &                                  &                            & Myanmar     & Southeast Asia \\ %\cline{2-6}
                       &                        &                          &                                  &                            & Philippines & Southeast Asia \\ %\cline{2-6}
                       \cline{1-7}
{\multirow{2}{*}{B422}}& {\multirow{2}{*}{27}}  & {\multirow{1}{*}{ST1}} & {\multirow{1}{*}{East Asian}}     & {\multirow{1}{*}{Beijing}}   & Malaysia       & Southeast Asia \\ %\cline{2-6}
                       &                        & {\multirow{1}{*}{ST38}} & {\multirow{1}{*}{Euro-American}} & {\multirow{1}{*}{X2}}        & Philippines    & Americas \\ %\cline{2-6}
                       \cline{1-7}
{\multirow{1}{*}{B425}}& {\multirow{1}{*}{2}}   & ST93 & Euro-American & LAM5 & Honduras & Americas \\ %\cline{2-6}
                       \cline{1-7}
{\multirow{6}{*}{B525}}& {\multirow{6}{*}{11}}  & ST167 & Euro-American & T1 & {\multirow{6}{*}{Haiti}} & {\multirow{6}{*}{Americas}} \\ %\cline{2-6}
                       &                        & ST42  & Euro-American & LAM9 & & \\
                       &                        & ST57  & Euro-American & LAM10-CAM & & \\
                       &                        & ST904 & Euro-American & T5 & & \\
                       &                        & ST904 & \emph{M. africanum} & AFRI\_1 & & \\
                       &                        & ST904 & \emph{M. africanum} & AFRI\_1 & & \\
                       \cline{1-7}
\end{tabular}
\caption{Biclustering results on data configuration 3,4,5 using CMTF\_PARAFAC\_ALS and CMTF\_Tucker\_ALS models via UBF on the coupled matrix-tensor. Bicluster B411 and B412 suggests that Euro-American strain ST908 and ST904 infects patients from four spatially close countries in Americas respectively. Bicluster B421 suggests that transmissive Beijing strain ST1 is wide-spread and infects patients from three different TB continents. Bicluster B422 groups patients from two neighbour countries, Malaysia and Philippines, who are infected with Beijing strain ST1 and X2 strain ST38.}
\label{BiclusteringResults_3_4_5}
\end{table}

If we use spatial distance matrix only, factorization via CMTF\_PARAFAC\_ALS results in 21 density-invariant biclusters, and we picked 5 most stable biclusters among them. The stability values of these biclusters range from 0.26 to 0.32. Table \ref{BiclusteringResults_3_4_5} shows 3 of these biclusters. Bicluster B411 suggests that Euro-American LAM2 strain ST908 infects patients from Dominican Republic, Puerto Rico, Trinidad Tobago, and Unites States, all from Americas. Notice how geographically close countries are collected together in a bicluster in the host-pathogen association analysis by incorporating spatial proximity into domain knowledge. Bicluster B412 suggests that Euro-American T5 strain ST904 infects patients from Ecuador, Haiti, Trinidad Tobago, and United States, which are again all in Americas. Bicluster B414 includes strains of both Bicluster B411 and B412, and combines the two common countries in these biclusters. It suggests that Euro-American T5 strain ST904 and Euro-American LAM2 strain ST908 infect patients from Trinidad Tobago and United States. When we factorize the coupled matrix-tensor via CMTF\_Tucker\_ALS in UBF, 6 density-invariant biclusters are found, and we picked 5 most stable biclusters among them, with average best-match score ranging from 0.09 to 0.50. Table \ref{BiclusteringResults_3_4_5} shows 3 of these biclusters. Bicluster B421 points out that transmissive ST1 Beijing strain is wide-spread, and it infects patients from Taiwan, Barbados, Dominica, Malaysia, Myanmar, and Philippines, which cover 3 different TB continents: East Asia, Americas, and Southeast Asia. This shows that, even if we use spatial proximity matrix to narrow down transmission events, transmissive ST1 strain is still associated with patients from multiple TB continents. Bicluster B422 contains 27 patients from Philippines and Malaysia, both from Southeast Asia, which are infected with ST1 and ST38 strains. Notice how these countries are grouped together using the spatial proximity matrix. Bicluster B425 consists of 2 patients from Honduras, both infected with Euro-American LAM5 strain ST93.

If we concatenate both genetic and spatial proximity distance matrices, factorization via CMTF\_PARAFAC\_ALS does not assign any density-invariant biclusters. When the coupled matrix-tensor is decomposed via CMTF\_Tucker\_ALS, we find 5 density-invariant biclusters, with average best-match score values ranging from 0.11 to 0.30. Table \ref{BiclusteringResults_3_4_5} shows one of these biclusters. Bicluster B525 contains 11 patients from Haiti which are infected with Euro-American strains ST167, ST42, ST57, ST904, and \emph{M. africanum} AFRI\_1 strains ST187 and ST1867. The full list of biclusters can be found in supplementary material. Notice that there is no order in stability of biclusters found using CMTF\_PARAFAC\_ALS and CMTF\_Tucker\_ALS. However, biclusters found using CMTF\_Tucker\_ALS are more biologically coherent. This shows that that high stability does not imply biological relevance.

\subsection{Incorporating distance, but not time}
Finally, in the last data configuration, we use genetic distance, spatial distance, but not time. This reduces the mutation path length and transmission path length, which increases the likelihood of mutation between the set of strains and transmission between the set of patients. To do so, we contract and sum the host-pathogen tensor along the time mode, and concatenate genetic proximity matrix and spatial proximity matrix. We bicluster spoligotypes and countries using CMMF\_ALS on this dataset in UBF. There are 17 density-invariant biclusters, and we picked the ones with average best-match score of 90\% and above. Full list of these biclusters are in the supplementary material. Table \ref{BiclusteringResults_6} shows 2 of these biclusters, B64 and B66. Bicluster B64 contains 3 patients from Bangladesh infected with Indo-Oceanic EAI5 strain ST1391 and Euro-American T1 strain ST447. Notice that EAI5 is a generic sublineage of Indo-Oceanic lineage, and T1 is a generic sublineage of Euro-American lineage, and they are both not-well-defined. Bicluster B66 contains 19 patients from Haiti infected with Euro-American, East Asian and \emph{M. africanum} strains. Haiti is an island next to Dominican Republic and immigrants of Haiti must have brought strains belonging to various lineages.

%suggests a strong association between Chinese patients and Indo-Oceanic EAI5 strain ST89. This association does not exist in the literature and it is a new finding, but it may not be biologically relevant while being a strongly connected stable bicluster. Notice that in these biclusters, genetically close spoligotypes are associated with spatially close patients, although none of the biclusters include more than one country of birth.

\begin{table}[!]
\centering
\begin{tabular}{|c|c||c|c|c||c|c|}
\cline{1-7}
{\multirow{2}{*}{\textbf{Bicluster}}} & {\multirow{1}{*}{\textbf{Number}}} & \multicolumn{3}{|c||}{\textbf{Spoligotypes}} & \multicolumn{2}{|c|}{\textbf{Countries}} \\ \cline{3-7}
                                      & {\multirow{1}{*}{\textbf{of patients}}} & \textbf{SIT no} & \textbf{Major lineage} & \textbf{Sublineage} & \textbf{Name} & \textbf{TB continent} \\ \cline{1-7}
{\multirow{2}{*}{B64}} & {\multirow{2}{*}{3}}    & ST1391 & Indo-Oceanic  & EAI5 & {\multirow{2}{*}{Bangladesh}} & {\multirow{2}{*}{Indian Subcontinent}} \\ %\cline{2-6}
                       &                         & ST58   & Euro-American & T1   & & \\
                       \cline{1-7}
{\multirow{11}{*}{B66}} & {\multirow{11}{*}{19}}  & ST1162 & East Asian    & Beijing & {\multirow{11}{*}{Haiti}} & {\multirow{11}{*}{Americas}} \\ %\cline{2-6}
                        &                         & ST168  & Euro-American & H3      & & \\
                        &                         & ST398  & Euro-American & LAM4    & & \\
                        &                         & ST57   & Euro-American & LAM10-CAM      & & \\
                        &                         & ST874  & Euro-American & S       & & \\
                        &                         & UST256 & Euro-American & H1      & & \\
                        &                         & ST541  & East Asian    & Beijing & & \\
                        &                         & ST1867 & \emph{M. africanum} & AFRI\_1      & & \\
                        &                         & ST822  & Euro-American & LAM9      & & \\
                        &                         & ST546  & Euro-American & X3      & & \\
                        &                         & ST3    & Euro-American & LAM2      & & \\
                       \cline{1-7}
\end{tabular}
\caption{Biclustering results on data configuration 6 using CMMF\_ALS via UBF on the coupled matrix-matrix. Bicluster B64 groups patients from Bangladesh who are infected with two strains of not-well-defined sublineages: Indo-Oceanic EAI5 strain ST1391 and Euro-American T1 strain ST58.}
\label{BiclusteringResults_6}
\end{table}

%\subsection{Biclustering for each major lineage}
\subsection{Host-pathogen association within each major lineage}
The six phylogeographic major lineages determined by CBN are established. Therefore, we subdivide the patient dataset based on six major lineages, and run UBF on each of them. We used data configuration 6, since it resulted in both stable and biologically relevant biclusters in the complete patient dataset. We found the most stable host-pathogen associations for each major lineage and reported their biological relevance.

\begin{table}[!]
\centering
\begin{tabular}{|c|c||c|c|c||c|c|}
\cline{1-7}
{\multirow{2}{*}{\textbf{Bicluster}}} & {\multirow{1}{*}{\textbf{Number}}} & \multicolumn{3}{|c||}{\textbf{Spoligotypes}} & \multicolumn{2}{|c|}{\textbf{Countries}} \\ \cline{3-7}
                                      & {\multirow{1}{*}{\textbf{of patients}}} & \textbf{SIT no} & \textbf{Major lineage} & \textbf{Sublineage} & \textbf{Name} & \textbf{TB continent} \\ \cline{1-7}
{\multirow{3}{*}{B712}}& {\multirow{3}{*}{5}}    & UST251 & Euro-American & S & {\multirow{3}{*}{Mexico}} & {\multirow{3}{*}{Americas}} \\ %\cline{2-6}
                       &                         & ST478  & Euro-American & X2   & & \\
                       &                         & ST1154 & Euro-American & LAM9  & & \\
                       \cline{1-7}
{\multirow{6}{*}{B732}} & {\multirow{6}{*}{9}}    & ST471  & East-African Indian  & CAS1-Delhi & {\multirow{6}{*}{China}} & {\multirow{6}{*}{East Asia}} \\ %\cline{2-6}
                        &                         & ST25   & East-African Indian  & CAS1-Delhi & & \\
                        &                         & ST381  & East-African Indian  & CAS1-Delhi & & \\
                        &                         & ST21   & East-African Indian  & CAS & & \\
                        &                         & ST203  & East-African Indian  & CAS & & \\
                        &                         & UST167 & East-African Indian  & EAI5 & & \\
                       \cline{1-7}
{\multirow{4}{*}{B733}} & {\multirow{4}{*}{11}}   & ST381  & East-African Indian  & CAS1-Delhi & {\multirow{1}{*}{China}} & {\multirow{1}{*}{East Asia}} \\ %\cline{2-6}
                        &                         & ST25   & East-African Indian  & CAS1-Delhi & {\multirow{1}{*}{Dominican Republic}} & {\multirow{1}{*}{Americas}} \\ %\cline{2-6}
                        &                         & ST21   & East-African Indian  & CAS & & \\
                        &                         & UST167 & East-African Indian  & EAI5 & & \\
                       \cline{1-7}
{\multirow{4}{*}{B741}} & {\multirow{4}{*}{222}}  & UST1   & East Asian  & Beijing & {\multirow{4}{*}{United States}} & {\multirow{4}{*}{Americas}} \\
                        &                         & ST190  & East Asian  & Beijing & & \\
                        &                         & ST1    & East Asian  & Beijing & & \\
                        &                         & ST2    & East Asian  & H1      & & \\
                       \cline{1-7}
{\multirow{1}{*}{B742}} & {\multirow{1}{*}{3}}    & ST265  & East Asian  & Beijing & {\multirow{1}{*}{China}} & {\multirow{1}{*}{East Asia}} \\
                       \cline{1-7}
{\multirow{2}{*}{B743}} & {\multirow{2}{*}{4}}    & ST1162   & East Asian  & Beijing & {\multirow{2}{*}{Haiti}} & {\multirow{2}{*}{Americas}} \\
                        &                         & ST1168   & East Asian  & Beijing & & \\
                       \cline{1-7}
{\multirow{7}{*}{B751}} & {\multirow{7}{*}{17}}   & ST325  & \emph{M. africanum}  & AFRI\_1 & {\multirow{7}{*}{United States}} & {\multirow{7}{*}{Americas}} \\
                        &                         & ST326  & \emph{M. africanum}  & AFRI\_1 & & \\
                        &                         & ST187  & \emph{M. africanum}  & AFRI\_1 & & \\
                        &                         & ST181  & \emph{M. africanum}  & AFRI\_1 & & \\
                        &                         & ST319  & \emph{M. africanum}  & AFRI\_2 & & \\
                        &                         & ST331  & \emph{M. africanum}  & AFRI\_2 & & \\
                        &                         & UST229 & \emph{M. africanum}  & AFRI\_2 & & \\
                       \cline{1-7}
{\multirow{2}{*}{B761}} & {\multirow{2}{*}{3}}    & ST479  & \emph{M. bovis}  & BOV    & {\multirow{2}{*}{Dominican Republic}} & {\multirow{2}{*}{Americas}} \\
                        &                         & ST481  & \emph{M. bovis}  & BOV\_1 & & \\
                       \cline{1-7}
{\multirow{2}{*}{B762}} & {\multirow{2}{*}{9}}    & ST409  & \emph{M. bovis}  & BOV\_2 & {\multirow{2}{*}{United States}} & {\multirow{2}{*}{Americas}} \\
                        &                         & ST683  & \emph{M. bovis}  & BOV\_2 & & \\
                       \cline{1-7}
\end{tabular}
\caption{Biclustering results on data configuration 6 using CMMF\_ALS via UBF on the coupled matrix-matrix for each major lineage. Bicluster B712 suggests that Mexican patients are likely to be infected with UST251, ST478, and ST1154 strains, given that the pathogen is a Euro-American strain. Biclusters B741 groups 222 US patients and shows that US patients are very frequently infected with Beijing strains, including the transmissive ST1 strain. Biclusters B761 and B762 suggest that, given that MTBC is an \emph{M. bovis} strain, it is more likely to infect a patient from Dominican Republic, if it is a BOV or BOV\_1 strain, and more likely to infect a US patient if it is a BOV\_2 strain.}
\label{BiclusteringResults_7}
\end{table}

Table \ref{BiclusteringResults_7} shows some of the most stable and biologically relevant biclusters. The full list of biclusters can be found in the supplementary material. Bicluster 711 of Euro-American lineage in the list of supplementary material contains 613 US patients infected with 60 different strains. Bicluster 712 listed in Table \ref{BiclusteringResults_7} suggests a strong association with Mexican patients and pathogens with three Euro-American strains, S strain UST251, X2 strain ST478, and LAM9 strain ST1154. Notice that all strains belong to different sublineages of Euro-American lineage. The average best-match score of this bicluster is 0.82. Bicluster 721 of Indo-Oceanic lineage listed in the supplementary material suggests an association between 40 Chinese patients and 16 different Indo-Oceanic strains, belonging to various sublineages. The stability value of 0.9621 suggests that this is a strong host-pathogen association.

Bicluster B732 listed in Table \ref{BiclusteringResults_7} contains 9 Chinese patients infected with CAS1-Delhi, CAS and EAI5 strains of East-African Indian lineage. Similarly, bicluster B733 suggests that patients from China and Dominican Republic are likely to be infected with the following East-African Indian strains: CAS1-Delhi strains ST381 and ST25, CAS strain ST21, and EAI5 strain UST167. 222 patients in bicluster B741 suggests that US patients are infected very frequently with the following Beijing strains: UST1, ST190, and the transmissive ST1 strain. Similarly, Bicluster B742 and B743 suggests that Beijing strains infect patients from China (ST265) and Haiti (ST1162, ST1168) respectively. Bicluster 751 shows that US patients are infected with AFRI\_1 strains ST325, ST326, ST187, ST181 and AFRI\_2 strains ST319, ST331, UST229 of \emph{M. africanum} lineage. Bicluster B761 suggests that patients from Dominican Republic are likely to be infected with BOV strain ST479 and BOV\_1 strain ST481 belonging to \emph{M. bovis} lineage. On the other hand, bicluster B762 suggests that US patients are infected BOV\_2 strains ST409 and ST683 belonging to \emph{M. bovis} lineage. These two biclusters suggest that, given an \emph{M. bovis} strain, it is likely to infect a patient from Dominican Republic if it is a BOV or BOV\_1 strain, whereas it is more likely to have infected a US patient if the strain is a BOV\_2 strain.


\section{Discussion and Conclusion}
We developed the Unified Biclustering Framework (UBF) to find host-pathogen associations in tuberculosis patients. To our knowledge, this is the first study to restate host-pathogen association analysis as a biclustering problem. UBF is flexible in the sense that distance and time can be added into domain knowledge of data analysis via coupled matrix-matrix and matrix-tensor factorization. This enables genome-phenome data fusion in one unsupervised learning framework.

Each bicluster refers to a possible host-pathogen association. We found statistically significant biclusters, some of which represent well-known host-pathogen relationships and some of which reveal new associations. For instance, bicluster B16 shows the well-known association of patients from Philippines and EAI2-Manila strains. Similarly, bicluster B741 shows that many US patients are infected with Beijing strain ST1, a well-known initiator of many outbreaks in the US. On the other hand, we also found new patient-strain relationships via genome-phenome data fusion by adding genetic proximity, spatial proximity and time into domain knowledge. For instance, bicluster B422 groups patients from two neighbour countries, Malaysia and Philippines, who are infected with Beijing strain ST1 and X2 strain ST38. Biclusters B761 and B762 suggest that patients from Dominican Republic are infected with BOV and BOV\_1 strains of \emph{M. bovis} lineage, whereas US patients are infected with BOV\_2 strains of \emph{M. bovis} lineage. Note that although we picked statistically significant biclusters, statistical significance does not imply biological relevance \cite{StatSigBioRelevance}. However, these new stable biclusters lead to new host-pathogen associations.

Host-pathogen association analysis can be extended by adding new patient and strain attributes. As future work, we will add RFLP and MIRU, two biomarkers of MTBC, into this analysis. In addition, we will add other patient attributes such as age group, ethnicity, homelessness and other risk factors of TB. We will also speed up UBF using line search in ALS-based coupled factorization algorithms. This will enhance both the speed and accuracy of coupled factorizations, which will lead to more accurate host-pathogen associations.



%%\newpage
%\section*{Appendix}

\section*{Acknowledgments}
\noindent We thank Dr. Minoo Aminian for providing the KBBN sublineage labels of strains. This work was supported by NIH R01LM009731.


% IEEEtranS
\bibliographystyle{ieeetr} % ieeetr - if you want the order to be the order of citation
\bibliography{UBFbib}

\end{document}
